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ABSTRACT 


A boundary integral formulation for the analysis of 
circular plate bending under lateral loads is developed using 
Green’s functions. The formulation specifically applies to 
annular plates with arbitrary boundary conditions. The plate 
bending solution in the plastic range is determined using a 
numerical method of incremental loading. A computer program to 
perform the required calculations was developed and is 
presented. Results for three case studies are included and 
compared with results obtained by other methods. Plate 
behavior in the elastic range is in excellent agreement with 
other analytical solutions, and in the plastic range is in 
reasonable agreement with published results obtained uSing a 
finite element method. 
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CHAPTER 1 


INTRODUCTION 


The solution of plate bending problems in the plastic 
range has been approached using finite element, and more 
recently, boundary integral methods (also called boundary 
element methods). Unlike the finite element method, which 
discretizes the inside of the plate into a number of small 
elements, the boundary integral method discretizes only the 
plate boundary, thereby reducing the order of the problem by 


one. 


Many different approaches have been used to formulate 
plate bending problems using boundary integrals, with perhaps 
the most attractive proposed by Stern‘l> and BezineS2”. These 
authors, working independently, developed a direct boundary 
integal method using Green’s functions to solve the general 
plate bending problem in the elastic range. This work has been 


extended by Moshaiov and Vorus‘3> to include Plastic behavior. 


Symmetry is commonly used to reduce the size of the system 
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of equations needed to analyze a symmetric structure. In the 
case of axisymmetrically loaded circular plates the problem 
becomes one-dimensional. For the linear elastic case, it has a 
closed form analytic solution. However, some authors have 
instead treated this circular plate problem as a two- 
dimensional problem, perhaps for the sake of demonstration 
(Kamiya and Sawki<4>)., They have used symmetry to allow 
discretization of a portion of the boundary instead of the 


entire boundary. 


This thesis develops a general formulation for solving 
annular plate bending problems using boundary integrals that 
reduces the problem to essentially a one-dimensional problem by 
the use of a "ring" type Green’s function. Such a formulation 
is most useful, inspite of the closed form solution, as it 
allows the analysis of plates with different boundary 
conditions and loads with one algorithm. Moreover, it permits 
a solution in the plastic range, which is not possible with a 
conventional analytic approach. In addition, the simplicity of 
the one-dimensional formulation given here has a unique 


significance for the teaching of boundary element methods. 


Butterfield‘®> has developed a boundary element 
formulation for the one-dimensional elastic beam bending 
problem. Here, a Similar approach is taken for the annular 
plate. It is also extended to include non-linear behavior 


using an incremental load method as outlined by Moshaiov and 
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Vorus. Numerical results for three different annular plate 
configurations are presented and compared to results obtained 
using other methods. A discussion of these results follows, as 


well as suggestions for future work in this area. 
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CHAPTER 2 


DEVELOPMENT OF THE GOVERNING DIFFERENTIAL EQUATION 


This chapter reviews the derivation of the governing 
differential equation for a circular plate subjected to a 
symmetrically distributed lateral load. The derivation is 


based on that provided by Timoshenko‘®, using Kirchoff’s 


assumptions. It is, therefore, applicable only to the linear- 
elastic case. Treatment of non-linear behavior is addressed in 
Chapter 5. 


2.1 Moment—-Curvature Relationships 


The first step in developing the governing differential 
equation for a circular plate is to find a relationship between 
moment and curvature. Curvature in the radial direction (K,) 
and the tangential direction (Kg) for a symmetrically loaded 
circular plate is found using geometrical considerations. K,, 
which is the inverse of the radius of curvature in the radial 
direction (R,), can vary as a function of the distance (r) from 


the plate center. Examining the plate of Figure 2-1 at a small 
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radial distance dr away from r, the slope (%) can be expressed 
in terms of the radial distance r and the deflection (w) as 


follows 


dw 
$ = = a (2-1-1) 
dr 


The radial radius of curvature is 


dr 
R = a 
dé 


giving an expression for the radial curvature in terms of w 


i dé d w 
K = — = ———— = = (2-1-2) 
= R dr dr° 
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Undeflected 
position 


Figure 2-1. Deflection, slope, and curvature relationships 
in circular plate bending 


Due to symmetry, the tangential radius of curvature is 
constant at any given radius r, and for small deflections can 


be expressed as 


ss 


o | 


giving an expression for the tangential curvature 


17 





l ¢ 1 dw 


K = — = — = = 
R r r dr 
To relate moment to curvature, consider the small plate 
element illustrated in Figure 2-2. Define M, and Mg, the 
radial and tangential moments per unit length respectively, as 


follows: 


h/2 

M = | go 2 dz 

r r 
-h/2 
CASS 

n/ 2 

Mo = | T5 Zz dz 
So)! 


where h is the plate thickness. Assuming the plate is thin and 
hence experiences a two-dimensional stress state, Hooke’s law 
yields the following expressions relating stresses to strains 


i e e 
oe = — >. Bo as v Eo 
Clay, } 


(D211) 


E 


G = _—_—_—_—_——— | € : + vi€ = | 
6 Cees) 6 r 


where E is Young’s modulus, v is Poisson’s ratio, and ¢e,© and 
eég& are the elastic radial and tangential strains, respec- 
tively. In Chapter 5, which discusses non-linear behavior, a 
distinction will be made between the elastic strain and the 
total strain. In the plastic range, the total strain includes 
elastic and plastic components, in accordance with the 


following expressions: 
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where the superscript p refers to the plastic strain 


components. Since Chapter 2 treats only elastic behavior, the 


total strain is simply equal to the elastic strain. 







i ——-—-—-—- —--— 
nes: 


Figure 2-2. Differential circular plate element. 


Substituting Eq. (2-1-4) into Eq. (2-1-3) gives 
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M = ——__ é +ve Za 
r ee r 8 
=n 2 
(2-1-5) 
Rey2 
E 
M = oe | co ee Zaz 
8 @Qaw) 6 r 
=|) 


Finally, from geometrical considerations, strain and curvature 


are related by 





dw 
e = K Z = = Ae 
dr 
C2] Lo) 
1 dw 
E = K 2 5 -—-— 2z 
e © rdr 


which, when substituted into Eq. (2-1-5), give the desired 


moment-curvature relationship for a circular plate 








dw v dw 
M = - D 7 
dr rdr 
(Bala 
2 
1 dw dw 
M = - D —— + »v 
“ r dr dr’ 
where D, the flexural rigidity, is given by 
Eh. 
D = oo (2-1-8) 
12 (l-v ) 


90 





2.2 Equilibrium Equations 


For the moment-curvature relationships to be useful, an 
equilibrium equations in terms of moments and shears must be 
found. The governing differential equation is then developed 
by substituting the expressions relating moment to curvature 


into the equilibrium equation. 


To find an equilibrium equation, consider first the 
laterally loaded circular plate element illustrated in Figure 
2-3, where Q, is the radial shearing force per unit length and 
q is the distributed load (force per unit area). The couples 


acting along edges ab and cd due to radial bending moments are 
ab: Mr r dé 


dM 
M + fi dr 
r 


dr 


cd: 





je ae ele dé 








The shear forces acting along the edges are 


ab: Qr dé 


Oe -—— dr 
r 


cd: r+ dr ds 





dQ | 
dr 

















wee Y/N ve 


G 

k 
‘Si 
@ 
ag 

a 

a 

Figure 2-3. Shear forces and moments acting on differentia. 


plate element. 


Because of symmetry, there are no shear forces acting 
along edges ad and bc, leaving only the tangential bending 


moments 


These couples can be resolved into components acting 
perpendicular to the axis ro (which are equal and opposite and 
thus cancel), and components acting along the ro axis which 


have a total magnitude of 


Jd fo 





Mo dr dé 


Summing up all the couples and neglecting the small change in 
shear force across the element (i.e. dropping the dQ terms) 


gives the following equilibrium equation 


dM 
r 


dr 


Mr + 





dr | r+ dr dg - M . r dé 


Beg ar ae + Q. r dé dr = 0 


This equation is further simplified by ignoring higher order 


terms to give 


GC - M fie i = 0 (2-2-1) 





Also, from equilibrium in the z direction 


dQ. 
—_—— = q (2-2-2) 
air 


2.3 Differential Equation for Circular Plate Bending 


Having found the desired equilibrium equations, the moment 
curvature relationships of Section 2.1 are now used to develop 
the governing differential equation for the circular plate 
bending problem. Specifically, substituting Eq. (2-1-7) into 
Eq. (2-2-1) results in a third order differential equation to 


describe the response of a circular plate to lateral loading 
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d>w l dw 1 dw Q 
- + = 3 _ a ee = = (2-3-1) 
dr crar me Ndr D 








which can be rewritten to provide an expression for shear 


dow l dw 1 dw 
+ — : _-_ — (2-3-2) 
dr rdr je tel 








Differentiating Eq. (2-3-1) with respect to r and using Eq. (2- 
2-2) results in a fourth order differential equation in terms 
of r, w, D, andgq 


d’w arly l dw 1 dw q 
i : = a (DoE 


dr r fie aac dr’ r > dr D 











which can be expressed in a more concise form using the 


biharmonic operator as 


ve) 


w= = (2-3-4) 


where 


fo ee eee (2-3-5) 
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CHAPTER 3 


A BOUNDARY INTEGRAL METHOD 


This chapter describes a boundary integral method, and 
how it can be used to provide a general method of solution of 
circular plate bending problems. A one-dimensional beam 


bending problem is used to illustrate this method. 


3.1 General Formulation 

Solution of Eq. (2-3-3) analytically as a boundary value 
problem by applying known boundary conditions is possible 
(Timoshenko). However, this approach requires special 
attention for each plate configuration and load distribution. 
This difficulty is compounded when considering problems 
involving plastic behavior which should be solved numerically 


due to the non-linear behavior in the plastic range. 


Therefore, a general method is sought wherein the same 
procedure can be used to solve a variety of problems, thus 


lending the problems to computer solution techniques. One such 


method, adopted for this work to examine circular plate 








bending, is the boundary integral method. 


Stern and Bezine, working independently, obtained a 
general formulation for plate bending problems in terms of 
boundary integrals. These integrals involve displacements, 
slopes, bending moments, and shears on the plate boundary. 
Both authors used the following reciprocal identity, which can 


be derived using Green’s identity: 


D | Wo vw - ow VW dQ = 
N 
oo t Mgt - tg" - wee par (3-1-1) 
r 
where 

Q = plate domain 
r = plate boundary 
D = flexural rigidity 
V4 ~ biharmonic operator 
w,WqG —- deflections 
a.8C. = slopes 
M,Mg - bending moments per unit length 
Q,Qg - shear forces per unit length 


To arrive at the boundary integral equation, wq is selected as 


a Green’s function. 





3.2 Beam Bending Problem 


Butterfield has presented a reduced form of Eq. (3-1-1) 
for the one-dimensional case of beam bending, which is 
developed in the discussion that follows. As a first step, 
Butterfield multiplies both sides of the governing differential 
equation for a beam by a second function and integrates over 
the beam’s length, giving 


L i 
| EI V ww. dx = | Aw, dx (3-2-1) 
0 


G 
where A is the distributed load per unit length, I is the area 
moment of inertia of the beam cross section, and the operator 


V* is defined for the beam case as 





Integrating the left hand side of Eq. (3-2-1) by parts 


several times, Butterfield develops the equation 


Recognizing that for the second function 


4 
EI V Wo = AG 


the preceding equation can be rewritten as 
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| —Wo Q + $$ M - M, 6 + Q. w (3-2-2) 


To gain some physical insight into Eq. (3-2-2) it is 
useful to recall Betti’s reciprocal work theorem. This theorem 
states that for an elastic structure subjected to two 
independent causes (e.g. loads), the total work done by the 
first cause in moving through the displacements resulting from 
the second cause is equal to the total work done by the second 
cause in moving through the displacements produced by the first 
cause‘/>, The theorem, as it applies to our beam bending 
problem, can be better understood by examining the examples of 


Figure 3-l. 
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b. M; 9 ji $i Mj 


Figure 3-l. Examples of Betti’s reciprocal work theorem 


In Figure 3-la, the product of load Pj and displacement 
6ij is equal to the product of load let and displacement 6 ji- 
Similarly, in Figure 3-lb, the product of moment M; and slope 


#i:j 1s equal to the product of moment Mj and slope #33. 


This concept can be applied to understand Eq. (3-2-2), 
with the i cause the actual loading or applied moment condition 
for the system of interest and the j cause the loading or 
applied moment condition of the second function. In this 
equation, the A term represents the external load in the system 
of interest, which when multiplied by the deflection wag 
described by the second function, gives a reciprocal work term. 
Similarly, the #q term is a slope for the second function, 


which when multiplied by the moment M for the system of 
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interest yields another work tern. This concept applies to all 
the terms of the equation, which are combined according to 


Betti’s theorm of reciprocal work to yield the given equality. 


For Eq. (3-2-2) to be useful for solving the beam bending 
problem, Butterfield chooses as the second function a Green’s 
function which is based on a singular loading condition; that 
is to say, a point load. If the point load of unit magnitude 
is applied at either x=0 or x=L, the second intergral term on 
the left hand side of Eq. (3-2-2) takes on the value of the 
boundary condition of the beam of interest. Usually, some 
attention should be given to cases where the integrals become 


Singular. 


This yields two equations. However, in general, beam 
bending problems involve four unknown boundary conditions. For 
example, if simple supports are specified on both ends, the 
moments and deflections at the two ends are known to be zero, 
while the shears and slopes are not known. Therefore, two 


additional equations are required. 


To obtain two more equations, derivatives of the first two 


equations are taken. This gives an additional Green’s 
function. Boundary conditions can now be obtained by solving 
the four equations simultaneously. Knowing the boundary 


conditions, the beam problem can be solved for interior points 


using Eq. (3-2-2) and appropriate derivatives with the point 
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load positioned at the location of interest. 


The formulation for the beam bending problem outlined 
above is a direct approach based on Green’s identity. 
Butterfield also develops a boundary integral equation using 
the more intuitive indirect method. For the circular plate 
formulation that follows in Chapter 4, only the direct method 


is used. 





CHAPTER 4 


DEVELOPMENT OF THE FORMULATION FOR CIRCULAR PLATES 


This chapter presents a general method for solving 
circular plate bending problems in the elastic range using 
boundary integrals. The resulting integral equations are 
analagous to those developed by Butterfield for the one- 
dimensional beam problem, which are discussed in Chapter 3. 


Treatment of non-linear behavior is discussed in Chapter 5. 


4.1 Description of Plate Geometry 


Figure 4-1 depicts a symmetrically loaded annular plate, 
and includes notation that will be referred to throughout the 
development of the integral equation and the discussion of the 


selected Green’s functions. 


Important symbols are defined as follows: 





W,WG 7: deflections 


$,%¢ = slopes 

Mr »Moer —- radial bending moments per unit length 
Mg,Mgg —- tangential bending moments per unit length 
Qr,Qgr - radial shear forces per unit length 

Qg,Q8qg - tangential shear forces per unit length 

a = plate outer radius 

b = plate inner radius 

h = plate thickness 

r = radius of interest 

q — distributed load per unit area 





Figure : Symmetrically loaded annular | -:te. 
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4.2 Development of Integral Equation 


Though the laterally loaded circular plate bending problem 
is two-dimensional, symmetry reduces the problem to essentially 
a one-dimensional type problem. Therefore, it is expected that 
the integral equation developed by Butterfield for the beam 
case can be adapted to the case of the circular plate. This 


equation (Eq. (3-2-2)) is repeated below: 


L 


| A Wo - w AG dx = 
0 


To obtain some indication of what the circular plate 
integral equation will look like, a non-rigorous method is 
first used to adapt Eq. (3-2-2) to the circular plate case, 


followed by a more rigorous mathematical development. 


In the non-rigorous approach, a first step is to replace 
the integral on the left hand side with an area integral and 
the line load per unit length A with the load per unit area q. 
Next, adjustments are made to reciprocal work terms on the 
boundaries (right hand side of Eq. (3-2-2)). Specifically, the 
shear terms are replaced with the product of the radial shear 
per unit length and the total length (2rr), observing that for 
symmetrical loading, the tangential shear Qg is 0. Moment 
terms are replaced with radial moments, since there are no 


tangential moments on the boundaries. These moments must also 
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be multiplied by the length (2xr) to give total moment. The 


resulting expression is 


Zo 8 
| | Geen q¢ rdrd@= 
0 b 
a 
| 2er | “We o ti eo M - Ma. $ 6+ Vere Ww | (4-2-]) 
b 
which reduces to 
a 
| qd W, W do fe ale = 
b 
a 
| r —W oe 3 eo M = May? + Qa,¥ (4-2-2) 





So far, it has only been asserted that the integral 
equation for the circular plate case should have a form similar 
to Eq. (4-2-2). To obtain the exact integral equation, a 


direct method is used. 


As a starting point for the direct method, consider the 


following expression: 


ah alee 
G 
U = radar 
2 2 


dr dr 








This expression can be written slightly differently as V 
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dw dw l1 dw dw 
U = —_% = EP - | ae — aes ; | r dr 
dr dr b b dr r dr dr 
and 
dw d“wo - a dw 1 d“w dw 
V = a =s _ | — — : + rodr 
dr dr b b dr r dr dr 
Recalling Eq. (2-3-2), it is possible to rewrite these 
expressions in terms of shear Q, 
dw, dw 
U = — es = 
dr dr : 
a a 
dw. Q 1 dw. dw 
| ae Pedi = | aes ee dr 
b dr D b redr dr 
and 
dw d“w, 2 
V = — ; r = 
dr dr nS 
* dw Q, a 1 dw dw 
| =< ear - | -— —* as 
‘ dr D b le folie (ol) G 


Setting U equal to V 
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and canceling some terms gives 














dw dw | dw a 
a we ok = — -—- sr dr = 
adres ar b b dr D 
2 a a 
dw d Wo | dw Qa, 
— a - —_ —s— eect 
ar ar b b dr D 





The two integral terms above are integrated by parts, with 


the operator of Eq. (2-3-5) used to simplify the expressions 





= dw. Q Q cs 7 
Ger id r 4 
— — er dr = w. r- = Vww.,r dr 
ayo fs S 
b b b 
a a a 
ate Qer Qer 
—_ — r dr = wren _ Vw.~.wer dr 
dr D D G 
b b b 


which, when substituted into the previous equation, yield 











dw Eva 2 Q 2 a 
G re 4 7 
— rt —-— wi r-_— + Vww.r dr = 
a? dr’ Gp : 
ae b b b 
dw dw 2 Q = = 
G Gr 4 
— r - wre + Vw.~.wr dr 
Adc D : 
r b b b 


After some simplification, this expression becomes 


a 








4 4 
| Vwrer Wo -~ V Wo row dr = 
b 
2 2 = 
Q.. Qa, dw d Wo dw dw 
Wa > - we + — 5 -— — : 
D D ar-ar dr dr 
b 
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The last equation bears a close resemblance to Eq. (4-2- 
2), except for sign differences and the fact that terms 
involving the second derivative of deflection have yet to be 
replaced by terms involving moment. To accomplish the latter, 
the following expression is added to the third term of the 
right hand side of the preceding expression, and subtracted 


from the fourth term 


Through Eq. (2-1-7) this gives 








a 
v - Vv dr = 
| wor Wo Wo row r = 
b 
1 a 
— | c Wo o = Qar a? Mor - ie M | (4-2-3) 
D 
b 
which can be rewritten using Eq. (2-3-4) as 
a 
| | qW. - Wg W r dr = 
b 
a 
| iG | Wo Qs. - Ww Qap + ¢ Ma. - # M (4-2-4) 
b 








Remarkably, with the exception of a sign difference that is 
attributable to the difference in sign convention, this 
expression is identical to Eq. (4-2-2), which was asserted 


based on the integral equation developed by Butterfield. 
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4.3 Selection of a Green’s Function 

As is the case for the beam described in Chapter 3, Eq. 
(4-2-4) is useful only if wq is chosen as a Green’s function 
with certain properties. As will become evident later in this 
chapter, the Green’s function must describe a ring loaded plate 
to solve circular plate bending problems, just as a Green’s 
function describing a point loaded beam was used to solve the 


beam bending problen. 


To solve the beam problem, Butterfield choses a Green’s 
function describing an infinitely long beam with a point load. 
It is not essential that a Green’s function for an infinite 
beam be used. In fact, as one would expect from Betti’s 
reciprocal work theorem, a Green’s function describing a finite 
geometry will work as well, provided the function also 


describes the response to a point load. 


Applying Butterfield’s method to the circular plate case, 
a Green’s function representing the ring loaded plate shown in 
Figure 4-2 is adopted. The Green’s function and associated 
derivatives are included in Appendix A. It should be noted 
that the outer and inner radii of the Green’s function (c and 
d) need not coincide with the outer and inner radii of the 


actual plate (a and b). 
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Figure 4-2. Plate geometry for the selected Green’s 
function 


That a simple support is chosen for the outer plate radius is 
not important; a plate clamped at both the inner and outer 
radii would also have been suitable. What is important is that 
the loading condition for the Green’s function case be a ring 
load. To simplify calculations, a ring load of magnitude 1 is 


used. 


4.4 Solution of Integral Equation 


Having selected a Green’s function, the next step is to 
obtain a form of Eq. (4-2-4) suitable to solve the actual plate 
bending problem. Rearranging terms, Eq. (4-2-4) can be 


rewritten as. 
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aE | q Wr dr (4-4-1) 
b 


The ring load may be expressed in terms of the dirac delta 


munction 6&(r,r,) 


G Ig 6(r,r_) 


Substituting this expression into Eq. (4-4-1) and making use of 


the following property of a dirac delta function 





co 
| 6(r,r) wor dr = r w(r=r) 
—o 
gives 
a 
r Siar) = | r | —Wo o a 2a M = Mor $ + Qcr Ww | 

b 

a 
+: | q wor dr (4-4-2) 

b 


By choosing ro at radius a, and then at radius b, the 
deflection w in the left hand side of Eq. (4-4-2) is the 
deflection at the boundary, which is either a known or unknown 
boundary condition. Hence, two equations are available to 


solve for the four unknown boundary conditions. 
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To obtain two more equations, 
differentiated with respect to r,. 


Green’s function, 


again at 
slope of 
known or 


complete 


Eq. (4-4-2) is 
This yields the second 
which is evaluated with ro, at a and then 


b. The second Green’s function corresponds to the 


the plate, which at r=a and r=b is again either a 
an unknown boundary condition. There is thus now a 


set of four equations to solve for the four unknown 

















boundary conditions. 
Eq. (4-4-2) can be rewritten as 
a 
ie 
mor -r .) = | — “Wo Je a ae sii = Mor $ + Qcr W | 
r 
oO b 
a 
1 
a: r, | qwor dr (4-4-3) 
b 
Differentiating this expression with respect to ro gives 
a 
dw(r=r_) r dw dé dM dQ 
- ov _ = a G Q + G M Gr aan Gr 
dr r dr dr dr dr 
oO re) b 
1 ss dw wine 
+ — | q — r dr = a (4-4-4) 
rs dr re 
Eqs. (4-4-3) and (4-4-4), evaluated both for ro=a and ro=b, 
give rise to a system of four equations that when solved 
simultaneously, yield the four unknown boundary conditions. In 


matrix forn, 


these equations can be written as 
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aw(a) “Wa (a, a) $,(a,a) —-M,(a,a) Q,,.(a,a) Q (a) 
bw(b ) - “wa(a,b) $,(a,b) “Mg. (a,b) Q,,. (a,b) M (a) 
—~a$(a) = (EY, $.(a,a) =o BRED Qa. (aya) $(a) 
-~b4#(b) -w (a,b) $-(a,b) “My (a,b) Gp (2+) | w(a) 
~wWo(b, a) $,(b,a) —Mg ,.(b, a) Q,,.(b, a) Q(B) 
4 [rMa(PeB) Fg (DsB) Mg p(b»b) Qg. (bb) | [M, (b) 
-wa(b,a) $.(b,a) -M,.(b, a) Q..(b,a) $(b) 
Jnwg (bb) $-(b,b) -Mo_.(b,b) Q..(b,b) w(b) 
a 
|, q Wo(r,a) r dr 
a 
| q Wo(r,b) r dr 
b 
= i. (4-4-5) 
| q wa (r,a) rodr _ w(a) 
b 
2 1 
| qwa(r,b) r dr -  w(b) 
b 
where 
| }(s,t) = | \(r=sjr st) 
and 
poy. a) 
dr 


Using Gauss elimination, the preceding system of equations is 


solved for the four unknown boundary conditions. 


The final step is to calculate deflection, slope, moment 


and shear across the plate. Deflection and slope are found 
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using Eqs. (4-4-3) and (4-4-4) by substituting in the 


calculated boundary conditions and evaluating the expressions 


at the value of ro of interest. Moment and shear are 


determined using Eqs. 


(2-1-7) and (2-3-2), which require 


calculation of the second and third derivatives or w with 


respect to ro. Differentiating Eq. (4-4-4) results in the 


desired expressions 


2 
d w(r=r_) ; 2 - 
alr * r 
O O 
l 
+ — 
r 
and 
d°w(r=r_) 
_ es 
dr r 
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a 
2 2 2 
d ie d Mor d Qar 
r M rr ¢ Saas ey 
dr dr dr 
‘e) O O 
b 
2 dw(r=r_) 
r dr = - -— (4-4-6) 
r dr 
O O 
3 3 3 2 
d ¢ d M d Q 
G M _ Gr $+ Gr a 
dr - dr dr 
b 
2 
3 d w(r=r) 
r dr - = ; (4-4-7) 
r dr 
O O 


making a complete solution for the circular plate bending 


problem in the elastic range possible. 
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CHAPTER 5 


NON-LINEAR BEHAVIOR 


This chapter develops a boundary integral formulation for 
circular plate bending problems that includes non-linear 
material behavior. Specifically, the boundary integral 
equation developed in Chapter 4 for the elastic plate is 
modified to include plastic moment terms to account for the 
plastic behavior. An incremental load method similar to that 
used by Moshaiov and Vorus can then be applied to arrive at a 


solution to the plate bending problem once yielding has begun. 


5.1 Plastic Stress-Strain Relationships 


When considering plasticity, strain can be thought of as 
having two components: an elastic component and a plastic 
component. In cylindrical coordinates, which is applicable to 
the two-dimensional circular plate case, the total strain e¢e can 


be expressed in terms of these two components as 





(S=f=)) 


where the superscripts e and p refer to the elastic and plastic 


strain components, respectively. 


Eq. (2-1-4), which provides the relationship between 
stress and strain for a linearly elastic material, can thus be 


rewritten in terms of the total strain and plastic strain as 


follows: 
E 
ee diy [Perth be ed 
(5=1— 2») 
E 
To = Gosey | | ar = re Te v | 3. = e | 


5.2 Plastic Moment-Curvature Relationships 


To develop the moment-curvature relationships with 


plasticity, Eq. (5-1-2) is substituted into Eq. (2-1-3) to give 


lay 2 
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r 
Ma = 
where the radial 


are defined as 


Recalling the relationships for total strain of Eq. 
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and tangential plastic 
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z dz My 
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dz 


CUa=2=2Z) 


dz 


(2=1516) anid 


performing the integration over the plate thickness yield the 


moment-curvature relationship for a circular plate with 
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plasticity 








, - : dw v dw 3 
ro dr’ ‘ 5 ake ; "r 
(5-2-3) 
1 dw vdw : 
a ee ee 


5.3 Governing Differential Equation With Plasticity 


Chapter 2 develops the governing differential equation for 
a circular plate in the elastic range by substituting the 
elastic moment-curvature relationships into the elastic 
equilibrium equation. In this section, a similar approach is 
used to develop the governing differential equation with 
plasticity. The difference here is that the moment-curvature 


relationships (Eq. (5-2-3)) include plastic terms. 


Substituting Eq. (5-2-3) into the equilibrium equation for 


a@ circular plate (Eq. (2-3-1)) gives 











dw i d‘w l dw l = dM - 
- + — iad ar aan 1 M + —— ri- Mo 
dr r dr r dr r D dr 
Q 
= _ (5-3-1) 
D 


where the tilda is used to distinguish the fact that the shear 
includes plastic effects. Differentiating this expression once 


with respect to r gives the fourth order governing differential 
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equation for the plate with plasticity 














dw 2 dw L d?w l dw 
ce + _— : r= aa ; + a, cae = 
dr r dr ro sd rr iar 
q 1 | 2 am P a’m P 1 am_P 
icc, ee - ae iG Sa pee 0 
D D rear dr’ roar 


which can be rewritten using the operator defined in Eq. (2-3- 
5) as 
Hen eed d°m P 1 aM_.P 

r 8 


q 
PO el! Oe - = (5-3-2) 
D D r dr dr r dr 





This equation is identical to the equilibrium equation for the 
elastic case with the exception of the terms that include 

derivatives of the plastic moments. Note that these terms may 
be thought of as additional pseudo loads that must be combined 


with the actual load to account for the plastic behavior. 


5.4 Boundary Integral Formulation With Plasticity 


In the preceding chapter, the governing differential 
equation is used to develop a boundary integral formulation for 
a circular plate in the elastic case. This section presents a 
similar approach to develop a boundary integral formulation 


for the plastic case. 


As a starting point, Eqs. (2-3-4) and (5-3-2) are 


substituted into Eq. (4-2-3) to give 





a 











| Iqwr dr = | r [-w, Q.. + eo M = Moy $ + Qar Ww 
b b 
(5-4-1) 
: e iezeane © cate el dm ,P 
+ qw,rerdr ~- = 37 -— — We i cdr 
3 | d dr° Eeday € 
b b r dr 
This equation is not yet in a form that is usable for 
solving the circular plate bending problen. In the first 


place, there are several terms that include derivatives of the 
radial or tangential plastic moments, which are not generally 
known across the plate. Secondly, the shear and moment terms 
Q, and Mr evaluated at the boundaries are not the actual shears 
and moments, since they were developed based on elastic 


considerations only. 


To resolve these concerns, the integral term involving 


plastic moments in Eq. (5-4-1) is integrated by parts as 


follows: 
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dM? d*m P dm,” 
ae 2— + 4 5 —-— dr = 
dr dr dr 
b 
a 
P dM? Pp 
Mo ae + r— = My 
dr 
b 
a 
dw dm P 
- —_ Mola Spo! Co liege 
dr e dr : 
b 





a 
dm P d’m_ P dM? 
Wo 2— + r , - — dr = 
dr dr dir 
b 
a 
dm ? dw 
Wo MF cee meat = ae - r— eel 
dr dr 
b 
a ; : 
W Ww 
+ a Mat + r = MF dr 
dr dr 
b 


Using the above relationship and choosing a Green’s function 
that describes a ring loaded plate as is done for Eq. (4-4-2), 


Eq. (5-4-1) becomes 
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a 
a abst.) 7 | = | me = is Fe He a Mee ee ee i | | 
b 
a 
dM? dw 
s abe MP oe a Sag a gee IM 
dr dr - 
b 
a a , 
dw dw 
== qw,rdr —- = MP 4 & rm P dr 
G 3 2 r 
dr dr 
b b 


Recall that Q@, and Mr represent shear and moment for the 
elastic case. However, they are not the actual shear and 
moment for the plastic case. With plasticity, the actual 
shears and moments are combinations of elastic and plastic 
terms. From Eq. (5-3-1), the shear relationship with 


plasticity is 


and from Eq. (5-2-3), the moment relationship with plasticity 


is 


where, as with the case for shear, the tilda is used to 
distinguish the fact that the moment term includes plastic 
effects. When substituted into the preceding expression, these 
expressions give the desired boundary integral formulation for 


the circular plate bending problem 


ay 











C 
ger =r o) 7 | GE ve " se ae = Wa x: Sh Se > | 
r 
oO b 
a 
it 
+ —_— q Wa © dr 
ae 
b 
2 
1 dw d Wo 
ae ee me MP + —s'rMP | ar (5-4-2) 
r 
ae dr dr 
b 


By evaluating Eq. (5-4-2) with the ring load positioned 
both at ro=a and ro=b, two equations are obtained to solve for 
the unknown boundary conditions. Eq. (5-4-2) is differentiated 
with respect to ry, to obtain a second Green’s function. This, 
via Eq. (2-1-1), gives an expression for slope of the actual 
plate. The second Green’s function is evaluated with the ring 
load positioned at ro=a and ro=b, to provide the remaining two 


equations. 


Differentiating Eq. (5-4-2) with respect to ro gives 
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dr r dr dr : dr dr 
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a 
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We q Be ed: = = 
ro dr r 
b oO oO 
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4|| d & 
l dr dr’ 
me ou a Mae + a r M ae (5-4-3) 
Po dr dr 
oO 
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In matrix form, the system of equations to be solved for 
the four unknown boundary conditions for the non-linear case 


can be written using the notation of Eq. (4-4-5) as follows: 
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Observe that Eq. (5-4-4) is identical to Eq. (4-4-5) 
except for the additional moment terms in the last matrix of 
Eq. (5-4-4). These additional terms account for the plastic 


behavior. 


To find moments and shears across the plate, expressions 


for the second and third derivatives of deflection with respect 
































to ro are required. These expressions follow: 
2 2 2 2 2 a 
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CHAPTER 6 


NUMERICAL SOLUTION 


This chapter describes the numerical procedure that was 
used to solve the circular plate bending problem in both the 
elastic and plastic ranges. The incremental load approach used 
in the plastic range is first described. A brief description 
of the computer program that was developed based on the 
equations of Chapters 4 and 5 follows, including a discussion 


of some aspects of numerical methods used in this progran. 


6.1 Incremental Load Method 

Once plastic deformation starts to occur in the plate, the 
linear relation between stress and strain is not valid. As 
discussed by Mendelson<®>, strains in the plastic range are no 
longer uniquely determined by the stresses, and in fact depend 
on the loading history. Therefore, it is not possible to use 
the relationships developed in Chapter 5 to directly arrive at 
a plate bending solution given a loading condition outside the 
elastic range. Rather, an incremental approach must be used 


wherein the load is increased in small steps once yielding 
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occurs at any point in the plate, and the complete stress state 


for the plate is determined before load again is increased. 


The incremental load method used for this analysis is 
based on that adopted by Moshaiov and Vorus, and is frequently 
seen in finite element solutions of structural analysis 
problems. The principal steps of this method are outlined 


below: 


l. Using the plate parameters and loading conditions for 
the bending problem at hand, an elastic solution is obtained 


with the equations developed in Chapter 4. 


2. The load at which the onset of yielding will occur in 
the plate is calculated based on the stress state determined in 
Step 1. The selected yield criterion is the von Mises 
criterion, which for the symmetrically loaded circular plate 


case has the form: 


where og is the equivalent, or effective, stress which 
represents the von Mises yield surface. Yielding will occur 


when og equals or exceeds the uniaxial yield stress (Sy). 


Gr. Unknown boundary conditions at yielding are 


determined using Eq. (5-4-4), with plastic moments initially 
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set equal to zero. Total strains at predetermined integration 
points are calculated using Eqs. (2-1-6), (5-4-3) and (5-4-5). 
From Eq. (2-1-4), the stress state at the integration points at 
yielding is determined, and is stored to be later updated as 
load is increased. Also stored for later updating are the 


deflections, slopes, moments and shears across the plate. 


4. The load which caused yield is then increased by a 
small incremental amount. The increments of the total strains 
resulting from this incremental load are calculated using Eq. 
(2-1-6) in an incremental form, from which elastic strains can 


be calculated by means of Eq. (5-1-1). 


5. Using Eq. (2-1-4), the elastic stress increment 
corresponding to the applied incremental load is calculated. 
This incremental stress is added to the stress stored in Step 3 
to give the total stress at the end of the load increment. A 
check is made throughout the plate using the yield criterion of 


Step 2 to determine where yielding has occurred. 


6. In those plate regions where yielding has occurred, 
the plastic strain increment (that is, the plastic strain due 
to the incremental load) is calculated. Only materials that 
exhibit elastic-perfectly plastic behavior as depicted in 
Figure 6-1 have been examined. However, the method can be used 
with other material behavior, such as strain hardening. For 


materials exhibiting elastic-perfectly plastic behavior, the 
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plastic strain increment 6&p 1s the total incremental strain 


determined from Eq. (2-1-6). This is shown in Figure 6-1. 


7. The plastic moments are determined using Eq. (5-2-2). 


8. Steps 3 through 6 are repeated, with the plastic 


moments calculated in Step 7 used as an input to Eq. (5-4-4). 
Step 7 1s repeated and the plastic moments compared to the 
moments calculated previously in Step 7. Iterations are 
performed as necessary until these values converge, according 


to a predetermined convergence criteria. 


9. After convergence of plastic moments is achieved in 
Step 8, stresses, plastic strains, deflections, slopes, moments 


and shears throughout the plate are updated. 


10. Another increment of load is applied, and Steps 3 
through 9 above are repeated, using the plastic moment from the 


previous load step as an initial estimate for plastic moment. 
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Figure 6-1. Uniaxial stress-strain curve for eiastic- 


perfectly plastic material. 


6.2 Computer Program Description 


Computer programs written in the Fortran 77 language have 
been developed to solve the circular plate bending problem and 
are included as Appendix B. Three programs, supported by nine 


subroutines, are used. 


The programs INPUT and LOAD create data files which are 
used by the program MAIN to solve the bending problen. INPUT 
and LOAD prompt the user for information describing the problem 
to be solved, such as material properties, plate geometry, 
boundary conditions and plate loading conditions. In addition, 


certain adjustable parameters are input, which include step 








sizes for the load increments, number of integration 
increments, and the desired percentage for convergence of the 
plastic moment increments. The program MAIN determines the 


plate bending solution, both in the elastic and plastic ranges. 


A flow diagram outlining the major elements and overall 


logic path of this program is shown in Figure 6-2. 





FUNCTION BLOCK PROGRAM/SUBRDUTINE 










INPUT 
PEFINE PROBLEM LOAD 
START MAIN PROGRAM MAIN 
INTEG 
AMATR 
SOLVE ELASTIC BENDING PROBLEM GAUSS 
FOR GIVEN LOADING CONDITION SEE 
SOLVER 
FIND LOAD AT YIELDING AND SDLVE YIELD 
ELASTIC BENDING PROBLEM 
STORE RESULTS FOR LATER UPDATING UPDATE 
IN THE PLASTIC RANGE 
INTEG 
APPLY INCREMENTAL LOAD cule 
AND SOLVE BENDING PROBLEM pegs 
SOLVER 
CALCULATE PLASTIC MOMENTS PLAST 
AND COMPARE TD ASSUMED MOMENTS 
ts 
CONVERGENCE CONDITION 
SATISFIED? 
ADD INCREMENTAL STRESSES, STRAINS, panes 
ETC TO PREVIOUS TOTALS 
RESUL 





DISPLAY AND STORE RESULTS 


is 
ULTIMATE STRAIN 
EXCEEDED? 


YES 


Figure 6-2. Computer program flow diagram. 
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6.3 Numerical Methods 

Eqs. (5-2-1), (5-4-4), (5-4-5) and (5-4-6) require 
integrations to be performed across the plate. These 
integrations are accomplished by a trapozoidal rule numerical 
integration scheme. It was found that 10 increments across 
both the plate radius and the plate half thickness gave 


satisfactory results. 


To ensure convergence of plastic moments, a root mean 
square average of the the plastic moments at each station 
across the plate is calculated and compared to the root mean 
Square average from the previous iteration. An agreement of 
less than 0.1% was used for the analysis work for which results 


appear in Chapter 7. 


65 








CHAPTER 7 


RESULTS AND DISCUSSION 


This chapter presents results of analyses performed on 
three different plate configurations using the computer program 
described in the preceding chapter. Results in the elastic 
range for all three cases are compared with the analytic 
solution from Roark*<9>, and for one case in the plastic range 
where published results using another solution method are 
available. Descriptions of the case studies are given in 
Sections 7.1 through 7.3. A discussion of the results is given 


in Section 7.4. 


7.1 Simply Supported Elasto-Plastic Annular Plate 


The annular plate shown in Figure 7-1 is first examined. 
The plate is simply supported at both the inner and outer 
radii, and is subjected to a uniform lateral load. A material 
exhibiting elastic-perfectly plastic behavior as depicted in 
Figure 6-1 is used. The yield stress (Sy) is 16 ksi and the 
Young’s modulus is 10x10°% ksi. Material properties and plate 


dimensions were selected to allow comparison with results 








obtained using a finite element method by Armen et al<1l9>, 


Plate deflection at yielding is shown in Figure 7-2, which 
also includes results from the analytic solution given by Roark 
for comparison. Plate deflections for the present solution at 


several loads in the plastic range are shown in Figure 7-3, as 


well as results obtained by Armen et al using a finite element 


method of analysis. 


Figure 7-4 shows the plastic zones in the plate at the 
three load conditions depicted in Figure 7-3. Results obtained 


by Armen et al are also included. 


y= 16 KSI 
EF = 10x10° KSI 
Y= 0.24 
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LO a 





Figure 7-1. Simply supported annular plate problem 
description. 
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‘igure 7-4. Plastic zones in simply supported annular pliets. 


7.2 Clamped Elasto-Plastic Annular Plate 


The annular plate shown in Figure 7-5 is next examined. 
The plate is clamped at both the inner and outer radii, and is 
again subjected to a uniform lateral load. A material 
exhibiting the same properties as described in Section 7.1 is 
used. For the clamped annular plate case, results obtained 
using other methods could not be found for the plastic range. 
Consequently, material properties and plate dimensions were 
selected to allow a qualitative comparison with results for the 


simply supported plate of Section 7.1 


Plate deflections at the onset of yielding are shown in 


Figure 7-6, and at several loads in the plastic range in Figure 
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7-7. The deflection at yielding using the analytic solution 
given by Roark is shown in Figure 7-6 for comparison. Finally, 
Figure 7-8 shows the plastic zones in the plate at the three 


load conditions depicted in Figure 7-7. 
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Figure 7-6. Deflection of clamped annular plate at the 
onset of yielding 
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Figure 7-7. Deflection of clamped annuler plate in the 
plastic range 
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Figure 7-8. Plastic zones in clamped annular plate. 


7.3 Simply Supported/Guided Elasto-Plastic Annular Plate 


The third annular plate to be examined is shown in Figure 


7-8. The plate is simply supported at the 


outer radius, with 


the inner radius guided. The plate is again subject to a 


uniform lateral load. A material exhibiting the same 


properties described in Section 7-1 is assumed. 


For this case, results obtained using 
again not be found for the plastic range. 
for a continuous simply supported circular 
obtained by Moshaiov and Vorus, and offers 


comparison with the results for the simply 


eZ 


other methods could 
However, a solution 
plate has been 
a qualitative 


supported/guided 


I INTR’DT OU <_—— 





annular plate. 


Plate deflection at the onset of yielding is shown in 
Figure 7-10 and at several loads in the plastic range in Figure 
7-1l for the present analysis method. The deflection at 
yielding obtained using the analytic solution from Roark is 
shown in Figure 7-10 for comparison. Figure 7-12 shows the 
plastic zones in the plate at the three load conditions 
depicted in Figure 7-11 for the present solution. Corres- 
ponding results from Moshaiov and Vorus for the continuous 


plate of Figure 7-13 are shown in Figures 7-14 and 7-15. 
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Figure 7-9. Simply supported/guided annular plat 
description. 
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Figure 7-10. Deflection of simply supported/guided annular 
plate at the onset of yielding. 
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Figure 7-11. Deflection of simply supported/guided annular 
plate in the plastic range 








7.4 Discussion of Results 

The results presented in Sections 7.1 through 1 wgaeeor the 
elastic range are uniformly in excellent agreement with the 
results using the analytic solutions given by Roark. This is 
as should be expected, since the elastic range solution using 


the present method is based on a closed form analytic solution. 


Results in the plastic range, while reasonable, are not in 
as close agreement with other published results. The only case 
where data for a direct comparison could be found is the 
configuration discussed in Section 7.1. From Figure 7-3, the 
deflections at loads of 652 psi and 852 psi for the present 
solution are in excellent agreement with the finite element 


results of Armen et al. 


Results for the annular plate of Section 7.1 at the 1052 
psi load are not as encouraging. Deflections and the extent of 
the plastic zone are greater for the solution obtained by Armen 
than for the present solution. To explain this difference, the 
size of the integration and load increments has been varied, 
but it shows little effect on the results. The reason behind 
this difference is not known. However, it is believed that it 
might be due to a programming error either in this report or in 
the reference. Also, differences could arise by not taking 
into account additional integration terms arising from the 


singularity at r=ro, as discussed in Appendix A. 








Finally, it is of interest to compare the results for the 
simply supported/guided plate of Figure 7-9 to the results for 
a simply-supported circular plate obtained by Moshaiov and 
Vorus. It is recognized that the two cases are not equivalent. 
However, one would expect the two plates to behave globally in 
a similar fashion. This is in fact the case, with both 


deflections and plastic zones showing reasonably good agreement 


in view of the differences between the two configurations. 








CHAPTER 8 


SUMMARY AND CONCLUSIONS 


8.1 Summary 


This thesis has developed a formulation for solving 
axXisymmetrically loaded annular plate bending problems using 
boundary integrals. This particular formulation is unique in 
that it treats the annular plate as a one-dimensional problen, 
using "ring" type Green’s functions to determine unknown 
boundary conditions and arrive at a plate bending solution. 
The formulation allows a numerical incremental load method to 
be used in the plastic range for the treatment of non-linear 
behavior. Results in the elastic range show excellent 
agreement with results obtained using a conventional analytic 
solution. In the plastic range, results are in reasonable 


agreement with those obtained using a finite element method. 


8.2 Conclusions 
This thesis treats the two-dimensional problem of analysis 
of axisymmetrically loaded annular plates as a one-dimensional 


problem. Using a method of solving one-dimensional problems 


ia 








with boundary integrals developed by Butterfield, the thesis 
demonstrates that a closed form solution for the annular plate 
problem can be obtained. It further demonstrates that this 
approach is suitable for the use of an incremental load method 


to obtain a solution in the plastic range. 


The formulation developed in this thesis is advantageous 
over a conventional analytic solution in that it allows a wide 
range of problems with different loading and boundary 
conditions to be solved using a single algorithm. Further, the 
simplicity of the approach is useful from an educational 


standpoint in the teaching of boundary element methods. 


8.3 Recommendations for Future Work 
There remains room for much additional work in this area. 


Some suggestions follow: 


l. A formulation for axisymmetrically loaded continuous 
circular plates remains to be developed. The work presented 
here lays a foundation for the continuous plate formulation. 
However, extending this method to the continuous plate case may 


require some modifications. 


2a In developing a formulation for the simple beam 
problem using boundary integrals, Butterfield uses both a 
direct and an indirect method, and shows that they yield the 


same result. Only the direct method is presented in this 
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thesis. The indirect method is a more intuitive approach than 
the direct method, and offers advantages in furthering the 
understanding of a boundary integral solution. It would, 
therefore, be worthwhile to pursue an indirect method as well 
as a direct method in developing a boundary integral 


formulation for circular plate geometries. 


3. The Green’s functions selected for this anlysis 
involve many terms and are awkward from a computational 
standpoint. More thought needs to be given to selecting a 
"ring" type Green’s function that is simpler and therefore 
offers advantages in simplifying calculations and reducing 


computational time. 


4. For the plastic range, the lack of close agreement 
between the results presented in this work and those obtained 
by Armen et al using a finite elemet method need to be better 
understood. Additional comparative data should be found or 
developed to increase confidence in the results presented in 


this work. 





APPENDIX A 


SELECTED GREEN’S FUNCTIONS 


This Appendix presents the selected Green’s functions and 
required derivatives that are used to solve the plate bending 
problems in this analysis. The functions which mathematically 
describe the deflection, slope, radial and tangential moments, 


and shear are obtained from Roark. 


A.l Description of Geometry and Sign Conventions 


As discussed in Chapter 4, the first Green’s function 
selected for this analysis describes the response of annular 
plate that is simply supported at the outer radius and 


unsupported along the inner radius to a ring load. Figure A-l 


depicts this plate configuration and important notation. 











Figure A-l. Representation of Green’s function used in 
the analysis. 


Roark uses a sign convention with deflection wa positive 
upward, slope ¢q positive when the deflection wg increases 
positively as r increases, moment M,; positive when creating 
compression on the top surface of the plate, and the shear 
force Qg;y positive when acting upward on the inner edge of an 
annular section. Subscripts c and d refer to the radial 


position. 


This sign convention differs from the sign convention of 
Timoshenko, which has been adopted for this analysis. 
Accordingly, adjustments must be made. Specifically, 
Timoshenko takes deflection to be positive downward and shear 


to be positive when acting downward on the inner edge of the 
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plate, such that the signs for these items as given by Roark 


must be changed. 


For slope, Roark uses a convention opposite to that of Eq. 
(2-2-1) and the convention used by Timoshenko. This 
difference, combined with the difference in defining the sign 
of deflection, effectively cancel each other, so the sign of 
slope from Roark does not change for this analysis. Finally, 
for moments, Roark and Timoshenko use the same sign convention, 


so no sign changes for moment terms are necessary. 


The formulas in the sections which follow include a number 
of general plate functions and constants that are not depicted 
in Figure A-l. Functions L and G and their derivatives are 
included in Section A.7; constants F and C are included in 


Section A.8. 


A.2 Singularities 


It is not possible to evaluate the Green’s function and 
associated derivatives at the exact radial location where the 
ring load is applied, due to singularities. To avoid this 
problem on the boundaries, the ring load is positioned a small 
distance inside the outer plate radius and outside the inner 
plate radius (0.00001 inches for the results presented in 
Chapter 7). During final review of this thesis, it was noted 
that because of the singularity at r=r,o, the domain integrand 


involving plastic moments of Eqs. (5-4-2) through (5-4-6) 








_ 


includes singularities. Therefore, the limiting values should 
be explored and a correction applied in way of the singular- 
ities. Insufficient time was available to repeat the analysis 
with the applied correction to see its effect on the results. 
Note that this correction affects only the solution in the 


plastic range. 


A.3 Deflection, Slope, Moment and Shear 


The following expressions describe deflection wa, slope 
$c, radial moment Mcy;y, and shear Qg corresponding to the first 


Green’s function: 





r 
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ee ou eden 7a cca: 2: 
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Weck 6) ol 3 








pc 
faa ae 
DC 
7 
D ~ plate flexural rigidity (Eq. (2-1-8)) 


Pg Magnitude of ring load 


0 
Si ea = singularity function 
The expression <r-r,>° represents a singularity function. 
[Mresitunction is equal to zero if r<r,. If r>ro, the brackets 
become like any other brackets. Hence, for r>rg, <r-rg?® is 


(r-ro) to the power of 0, which equals to l. 


A.4 First Derivative of Deflection, Slope, Moment and Shear 


The expressions which follow describe the first deriv- 
atives with respect to ro of wo, %, Mer and Qgc,r- The first 
derivative of wa corresponds to the second Green’s function 


required by the analysis. 
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A.5 Second Derivative of Deflection, Slope, Moment and Shear 


The following expressions describe the second derivatives 


with respect to ro of wa, %q, Mgr and Qgr 

















= 4 = Gd _ Be a F i.e 3 
dr dr dr ar 
oO oO 
2 2 2 2 
d ae d Fog r d G. 
= Z Fy i 0G i etre 
dr dr ar 
e) 
2 2 Z 
fn Gr = Geica 2 F —_ p r a 
dr’ dr’ 'G i “ dr ; 
oO oO 


87 








where 


A.6 


with respect to ro of wo, 
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Third Derivative of Deflection, Slope, Moment and Shear 


The following expressions describe the third derivatives 
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A.9 Derivatives of the Green’s Function With Respect to r 


The equations describing non-linear plate behavior 
included in Section 5-4 are in terms of the first and second 
derivatives of the Green’s function wq with respect to r. 

These derivatives can be expressed with the aid of Eqs. (2-1-1) 


and (2-1-7) in terms of moments and slopes as follows: 





dw 7 ; 
ane ic 
dr 
a’ M v 
ae 7" ear ig is 
dr* D r : 


Expressions for slopes and moments for the Green’s 
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function (¢q and Mg,) are available and have been presented in 
Section A.3. It is further noted that Eqs. (5-4-2) through (5- 
4-6) require derivatives of the above two expressions with 
respect to ro. This is accomplished via the expressions for 
derivatives of $qg and Mar with respect to ro included in 


Sections A.5 through A.7 above. 
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APPENDIX B 


COMPUTER PROGRAMS 


This Appendix contains the computer programs developed to 
perform the analysis work. A flow chart showing the overall 


program structure is presented in Chapter 6. 


The Appendix is organized with the programs INPUT and LOAD 
presented first, followed by the program MAIN. The supporting 


subroutines for the program MAIN follow that program. 








C 


Cx 

CEREKEEKESEKR EERE EKE TERETE AKER AKER KAKAKA TEAK TAREE AERA KEKEK EES 
Cx 

Cx PROGRAM INPUT 

Cx 

Cx THIS PROGRAM ALLOWS THE USER TO DEFINE THE PROBLEM OF INTEREST 


C* AND TO SET CERTAIN ADJUSTABLE PARAMETERS SUCH AS NUMBER OF DEPTH 
Cx INTEGRATION INCREMENTS. LOADING CONDITIONS ARE INPUT USING THE 

C* PROGRAM LOAD. TO RUN, THIS PROGRAM MUST ACCESS THE FILES GEOM.DAT, 
Cx BC.DAT,BCPOS.DAT AND ADJ.DAT 


Cx 
CHEEEEKEREEEREKAKARA TAKE EAE EAE AERA KATA E EERE REAR AKAAK AEE ERK EKE KS 
Cx 
C 

PROGRAM INPUT 
C 

IMPLICIT REAL*¥8 (A-H,0-Z) 

CHARACTER*1 BFLAG1(4),BFLAG2(4), RESP 

CHARACTER*19 LBL(8) ,LGEOM(9) 

CHARACTER*25 LADJ(6&) 
C 

REAL*8 ADJ(7),BC(4),GEOM(3),NU 
C 

INTEGER POS(8),BCPOS(8),M,N 
C 
Cx 
CEEEKEKEEKEKKAKEEEEKARE EEA ERE RE REE ERE KAKA AAA AEA AKA KEE ERE RAKE REESE 
Cx 
Cx INPUT MATERIAL CONSTANTS AND PLATE GEOMETRY. THE FLEXURAL 


Cx RIGIDITY IS CALCULATED BASED ON THE INPUT PARAMETERS. N IS THE 
C* TOTAL NUMBER OF PARAMETERS 


Cx 
OS SSS SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSsss sss sss sss sssssss: 
Cx 
C 
N=9 
C 
LGEOM(1)=’OUTER RADIUS A= ’ 
LGEOM(2)=’ INNER RADIUS B= ’ 
LGEOM(3)=’POISSON RATIO NU = ’ 
LGEOM(4)=’PLATE CONSTANT D = ’ 
LGEOM(5)=’*THICKNESS T = ’ 
LGEOM(6)=’YIELD STRESS = ’ 
LGEOM(7)=’YIELD STRAIN = ”’ 
LGEOM(8)=’ULTIMATE STRESS = ”’ 
LGEOM(9)=’ULTIMATE STRAIN = ’ 
C 
WRITE(*,1) 
] FORMAT(///,6X, ’WOULD YOU LIKE TO SEE THE CURRENT PLATE DIMENSIONS 
ts »/,2X,’ AND MATERIAL CONSTANTS (Y/N) 7? °\) 
READ(*,30) RESP 
C 


TPY (RESP | BO. 8° y¥” ) THEN 
WRITE(%,2) 
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2 FORMAT(//,2X,’THE CURRENT VALDES FOLLOW. NOTE THAT A ZERO’, 
+ /,2X,’ VALUE IS SHOWN FOR D WHEN T AND E ARE GIVEN AND VICE’, 
~ Juek, VERSE ..0°,/) 
OPEN(6,FILE=’GROM.DAT’,STATUS=’OLD’ ) 
WRITE(*,30) ° ° 
DO 8 I=1,N 
READ(6,70) GEOM(I) 
WRITE(*,6) LGEOM(I),GEOM(I) 


6 FORMAT(2X,A22,F20.10) 
8 CONTINUE 
CLOSE(6) 
ENDIF 
C 
WRITE(*,9) 
3 FORMAT(//,6X,’°DO YOU WISH TO INPUT OR CHANGE GEOMETRY OR’, 
1 /,2X,’MATERIAL CONSTANTS (Y/N) ? ’,\) 
READ(*,30) RESP 
C 
IF (RESP .EQ. ’Y’) THEN 
C 


WRITE (*, 10) 
10 FORMAT(/,2X, ‘INPUT THE VALUES USING DECIMAL NOTATION ....’,/) 
DO 12 I1=1,3 
WRITE(*,20) LGEOM(I) 
READ(*,70) GEOM(I) 
12 CONTINUE 


DO 13 1=5,9 
WRITE(4%,20) LGEOM(TI) 
READ(*,70) GEOM(I) 
13 CONTINUE 
GEOM(4)=GEOM(6) /GEOM(7) *GEOM(5)**3.0/12.0/(1.0-GEOM(3) **2.0) 


C 

OPEN(6, FILE=’GEOM. DAT’ , STATUS=’ NEW’ ) 

DO 15 I=1,N 

WRITE(6,70) GEOM(I) 
15 CONTINUE 
CLOSE(6) 
ENDIF 

C 
C 
SSS SSS SSS SSS SSS SSF FSS SSS ESE F FSFE FESS SSE FELL EES ESSE SSS SS SSS SS SoS 8 | 
C 
C INPUT PLATE BOUNDARY CONDITIONS 
C 
oF SSS SSS SSS SSS SSS SS SS SS SS SPSS FESS SESES SESE SEE E SELES SESS SESS SS SSS? Sd 8 8 | 
C 
C 


LBL(1)=’SHEAR AT A 
LBL(2)=*MOMENT AT A 
LBL(3)=’SLOPE AT A 
LBL(4)=’DEFLECTION AT A 
LBL(5)=’*SHEAR AT B 
LBL(6)=’MOMENT AT B 
LBL(7)=’SLOPE AT B 
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LBL(8)=’DEFLECTION AT B = ° 


OPEN(14,FILE=’BCLBL. DAT’ ,STATUS=’ NEW’ ) 
DO 155 1=1,8 
WRITE(14,156) LBL(I) 
155 CONTINUE 
156 FORMAT(A22) 
CLOSE( 14) 


WRITE(*,160) 

160 FORMAT(//,6X,’ WOULD YOU LIKE TO SEE THE CURRENT PLATE BOUNDARY’ 
+ »>/,2X,*°CONDITIONS (Y/N) ? ’\) 
READ(*,30) RESP 


IF (RESP .EQ. *Y’) THEN 
OPEN(11,FILE=’BC.DAT’,STATUS=’OLD’ ) 
OPEN(12,FILE=’BCPOS. DAT’ ,STATUS=’OLD’ ) 
WRITE(*,30) * ° 
DO 165 I=1,4 

READ(11,70) BC(T) 
165 CONTINUE 
DO 190 I=1,8 
READ(12,100) BCPOS(I) 
IF (BCPOS(I) .GT. 0) THEN 
WRITE(*,186) LBL(I),BC(BCPOS(I)) 


186 FORMAT (2X,A22,F20.10) 
ENDIF 
190 CONTINUE 
CLOSE(11) 
CLOSE(12) 
ENDIF 
C 


WRITE(*,17) 

17 FORMAT(//,6X,’°'DO YOU WISH TO INPUT OR CHANGE PLATE BOUNDARY’, 
+ Vier, CONDITIONS (Y/N) 2°, \) 
READ(*,30) RESP 


IF (RESP .EQ. ’Y’) THEN 
WRITE(*,19) 
HN) FORMAT(//,6X, INDICATE WHETHER THE BOUNDARY CONDITIONS FOR THE’, 
@ /,2X,°* PLATE ARE KNOWN OR UNKNOWN BY TYPING IN THE LETTER “"U"’, 
@ /,2X,°FOR UNKNOWN BOUNDARY CONDITIONS AND HITTING RETURN FOR’, 
@ /,2X,’ KNOWN BOUNDARY CONDITIONS’ ,//) 


WRITE(*,20) ’SHEAR AT A? ° 
READ(*,30) BFLAG1(1) 

WRITE(*,20) "MOMENT AT A? '° 
READ(¥,30) BFLAG1(2) 

WRITE(*,20) ’SLOPE ATA? ° 
READ(*,30) BFLAG1(3) 

WRITE(*,20) ’DEFLECTION AT A? °” 
READ(*,30) BFLAG1(4) 


WRITE(*,20) "SHEAR AT B ? °’ 
READ(*,30) BFLAG2(1) 


a 








C 


C 


20 
30 


40 


50 


60 


70 


WRITE(*,20) *MOMENT AT B 7? ° 
READ(#,30) BFLAG2(2) 

WRITE(4,20) “SLOPE AT B? * 
READ(%*,30) BFLAG2(3) 

WRITE(*,20) "DEFLECTION AT B ? ” 
READ(*,30) BFLAG2(4) 
FORMAT(2X,A22,\) 

FORMAT( Al) 


WRITE(*,40) 
FORMAT(//,6X, ’’TYPE IN THE KNOWN BOUNDARY CONDITIONS USING ’,/, 
@2X,** DECIMAL NOTATION WHEN PROMPTED’ ,//) 


OPEN(8, FILE=’BFLAG].DAT’ ,STATUS=’ NEW’ ) 
OPEN(9, FILE=’BFLAG2. DAT’ ,STATUS=’ NEW’ ) 
OPEN(10, FILE=’BC.DAT’,STATUS=’NEW’ ) 


DO 50 I=1,4 
WRITE(8,30) BFLAGI(1) 
IF (BFLAG](I) .NE. ’U’) THEN 
L=L+1 
WRITE(%,20) LBL(TI) 
READ(*,70) BC(L) 
WRITE(10,70) BC(L) 
ENDIF 
CONTINUE 


DO 60 I=1,4 
WRITE(9,30) BFLAG2(I) 
IF (BFLAG2(I) .NE. ’U’) THEN 
L=L+] 
WRITE(%,20) LBL(I+4) 
READ(x,70) BC(L) 
WRITE(10,70) BC(L) 
ENDIF 
CONTINUE 


FORMAT(F20.10) 


CLOSE(8) 
CLOSE(9) 
CLOSE(10) 


OPEN(11,FILE=’POS.DAT’ ,STATUS=’ NEW’ ) 
OPEN(12,FILE=’BCPOS.DAT’,STATUS=’ NEW’ ) 


L=0 
DO 80 I=1,4 
TF (CBELAGE( 1)... FQ.5°U"*) THEN 
L=L+] 
POS(I)=L 
BCPOS(I1)=0 
ELSE 
M=M+] 
BCPOS(I1)=M 
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80 


90 


100 


c 
Cx 
Cx** 
Cx 
Cx 
Ct 
Cx 
Cx 
Cx 
Cx 
Cx 
Cx 


POS(I)=0 
ENDIF 


WRITE(11,100) POS(TI) 
WRITE(12,100) BCPOS(TI) 


CONTINUE 


DO 90 I=5,8 

IF (BFLAG2(I-4) .EQ. °U’) THEN 
L=L+1 
POS(I)=L 
BCPOS(1)=0 

ELSE 
M=M+] 
BCPOS(I)=M 
POS(1)=0 

ENDIF 


WRITE(11,100) POS(I) 
WRITE(12,100) BCPOS(I) 


CONTINUE 
CLOSE(11) 
CLOSE(12) 
FORMAT( 11) 
ENDIF 


PE SETTT TTT TTT TTT TTT TEST TTT TTS SSTCSSSSCSSSSSSCSCTSSSSSCS SSS CTT TT 


INPUT CALCULATION ADJUSTMENT ITEMS. N IS THE NUMBER OF 
INPUT PARAMETERS IN THIS SUBMODULE. OBSERVE THAT SUBSEQUENT TO 
THE DEVELOPMENT OF THIS MODULE, PARAMETER NUMBER 3 WAS NO LONGER 
ACCESSED FROM ANY OTHER PROGRAM, AND THEREFORE SERVES ONLY AS A 
DUMMY PARAMETER. THE REMAINING PARAMETERS HAVE THE FOLLOWING 
MEANINGS 


1. EPSILON REFERS TO THE DISTANCE BETWEREN THE RING LOAD 
AND THE BOUNDARY. NOTE THAT IN SOME INSTANCES, LETTING 
EPSILON RQUAL TO ZERO CAUSES A DIVIDE BY ZERO TYPE ERROR. 


2. THE INCREMENTS FOR R REFER TO THE STATIONS ACROSS 
THE PLATE. THIS APPLIES TO THE NUMBER OF INTEGRATION 
INCREMENTS FOR DISTRIBUTED LOADS, AS WELL AS THE NUMBER 
OF DATA POINTS FOR THE RESULTS. 


3. THIS PARAMETER IS NO LONGER USED 


4. THIS REFERS TO THE EXTENT BY WHICH THE GREENS 
FUNCTION PLATE GEOMETRY OVERHANGS THE ACTUAL PLATE. 
FOR EXAMPLE, A VALUE OF 0.5 MEANS THAT THE INNER RADIUS OF 
THE RING LOADED PLATE DESCRIBED BY THE GREEN’S FUNCTION 
OVERHANGS THE ACTUAL PLATE BY 0.5 UNITS. 
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Cxre 
cz 
C 


102 


106 
108 


120 


as 


>. THE LOAD INCREMENT IS THE PERCENTAGE OF THE 
ORIGINALLY SPECIFIED LOAD MAGNITUDE THAT IS APPLIED DURING 
EACH LOAD INCREMENT AS THE PLATE IS LOADED IN THE PLASTIC 
RANGE 


6. THIS IS THE NUMBER OF INTEGRATION INCREMENTS USED TO 
DETERMINE THE PLASTIC MOMENT ACROSS THE PLATE. 


7. RMS IS THE ROOT MEAN SQUARE OF THE PLASTIC MOMENT 
INCREMENTS FOR EVERY POINT IN THE PLATE WHERE THE PLATE IS 
NO LONGER ELASTIC. THIS IS COMPARED WITH THE RMS VALUE FOR 
THE PREVIOUS INCREMENT, TO SEE WHETHER THE SOLUTION HAS 
SUFFICIENTLY CONVERGED. 1% CONVERGENCE MEANS THESE NUMBERS 
AGREE TO WITHIN LESS THAN 1% 


2 SPSS P SSE SSSSFLS SSS SSSSSCP ISP LSPISSS SOLIS SSS SCSLLIPES SESS SEES SSDS 2 FSS @ F © © ® | 


N=7 

LADJ(1)=’EPSILON FOR RO = ’ 
LADJ(2)=’INCREMENTS FOR R = ’ 
LADJ(3)=’WIDTH INTG INCR = ’ 
LADJ(4)=’BOUNDS FOR A AND B = ’ 
LADJ(5)=’LOAD INCREMENT = ’ 
LADJ(6)=’DEPTH INTG INCR = ’ 
LADJ(7)=’% RMS CONVERGENCE = ’ 


WRITE(*,102) 
FORMAT(//,6X, "WOULD YOU LIKE TO SEE THE CURRENT ADJUSTABLE’ 
+ ,/,2X,° PARAMETER SETTINGS? ’\) 
READ(*,30) RESP 
IF (RESP .EQ. ’Y’) THEN 
OPEN(16,FILE=’ADJ.DAT’,STATUS=’OLD’ ) 
WRITE(*,30) ’ ’ 
DO 108 I=1,N 
READ(16,70) ADJ(I) 
WRITE(*,106) LADJ(1I),ADJ(1) 
FORMAT(2X,A25,F20. 10) 
CONTINUE 
CLOSE(16) 
ENDIF 


WRITE(*,110) 

FORMAT(//,65X,’DO YOU WISH TO INPUT OR CHANGE ANY COMPUTATIONAL’, 
a /,2X%,° ADJUSTMENT ITEMS (Y/N) ? ’°,\) 

READ(*,30) RESP 


IF (RESP .EQ. ’Y’) THEN 


WRITE(%,120) 
FORMAT(/,’ INPUT THE VALUES USING DECIMAL NOTATION ....’,/) 
DOviZs f=1,N 
WRITE(*#,126) LADJ(1) 
READ(*,70) ADJ(1) 
CONTINUE 
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126 FORMAT(2X,A25,\) 
OPEN(13,FILE=’ADJ.DAT’,STATUS='NEW’ ) 
DO 130 I=1,N 
WRITE(13,70) ADJ(1) 
130 CONTINUE 
CLOSE(13) 
ENDIF 


WRITE(*,135) 
135 FORMAT(///,’ YOU ARE LEAVING THE INPUT MODULE’,///) 


END 
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C 


8 FSPPESESSSSSSSSSESSSSSISSS SS SISSIES SSSSSSSSSSSSSSSSSSESSSSSSSSSSS ESS SSS SS: 


Cx 
cx 
Cx 
Cx 
Cx 
Cx 
Cr 
Cx 
Ct 
Cx 
Cx 
Cx 
Cx 


PROGRAM LOAD 


THIS PROGRAM ALLOWS THE USER TO INPUT THE DESIRED LOADING 
CONDITION FOR THE PLATE. UNIFORM DISTRIBUTED LOADS, RAMPS, AND 
PARABOLICALLY DISTRIBUTED LOADS ARE PERMITTED, AS WELL AS RING 
LOADS. TO RUN, THIS PROGRAM MUST ACCESS THE FILE LOAD. DAT. 


USING THE INPUT DATA, THE PROGRAM CALCULATES THE TOTAL LOAD 
PER UNIT LENGTH THAT ACTS ON A RING OF THE PLATE SURFACE OR A 
WIDTH DETERMINED BY ADJUSTABLE PARAMETER #3 OF THE INPUT 
PROGRAM 


CEEETETERELEKE RTA EKSTRA ET ERE REREEEEAR EK ERE TERE ERE KARA EEA EKER 


Cx 
c 


C 


10 


is) 


30 
40 


PROGRAM LOAD 


IMPLICIT REAL*8 (A-H,0-Z) 
CHARACTER#1 RESP 
CHARACTER*30 LBL(4) 
REAL*¥8 LD(50) 

INTEGER N 


LBL(1)=’MAGNITUDE AT OUTER RADIUS = ’ 
LBL(2)=’DISTANCE FROM PLATE CENTER 


~ LBL(3)=’MAGNITUDE AT INNER RADIUS = ° 


LBL(4)=’DISTANCE FROM PLATE CENTER = ° 


OPEN(6,FILE=’LOAD. DAT’ ,STATUS=’OLD’ ) 
READ(6,90) LD 

N=INT(6.0+2.0*%LD(6) ) 

CLOSE (6) 


WRITE( 24,10) 
FORMAT(//,’ DO YOU WISH TO SEE THE LOADING CONDITIONS’, 
+ fi me ND? 7) 
READ(*,60) RESP 
IF (RESP .EQ. °Y’) THEN 
WRITEGE, 15) 
FORMAT(//,’ THE EXISTING LOADING CONDITIONS FOLLOW ...’,//) 
DO 40 I=1,4 
WRITE( 2,30) LBL(I),LD(I) 
FORMAT(2X,A30,F20.10) 
CONTINUE 
WRITE(%,60) ° ’ 
IF (LD(5) .EQ. 1.0) THEN 
WRITE(%,80; *LOADING IS LINEAR’ 
ELSE 
WRITE(#,80) ’'LOADING IS PARABOLIC’ 
ENDIF 
WRITE(*,60) ° ’ 
J=5 
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DOPe4oul=-7,6+INT(LD(G) ) 


JzJ+2 
WRITE(*,60) ” ’ 
WRITE(*,43) ’MAGNITUDE OF RING LOAD’,I-6, ”’ = *°,LD(J) 
WRITE(%,43) "’LOCATION OF RING LOAD’,I-6, * = °,LD(J+1) 
43 FORMAT(2X,A24,13,A43,F20.10) 
45 CONTINUE 


ENDIF 


WRITE(*, 46) 

46  FORMAT(//,2X,’DO YOU WISH TO CHANGE THE LOADING CONDITIONS’, 
Sue ec sN) 2 7X) 
READ(*,60} RESP 


IF (RESP .EQ. ’Y’) THEN 
WRITE( 2,50) 
50 FORMAT(//,2X,’DO YOU WISH TO CHANGE THE DISTRIBUTED LOAD’, 
oar CVI Ns) 
READ(*,60) RESP 
60 FORMAT(A1) 


IF (RESP .EQ. ’Y’) THEN 


WRITE(%,70) 
70 FORMAT(//,2X, SPECIFY THE LOAD PROFILE USING DECIMAL’, 
+ ie NOTATION cise 7? 5/7) 
DO 75 I=1,4 
WRITE(*,80) LBL(I) : 
READ(*,90) LD(I) 
75 CONTINUE 
80 FORMAT (2X, A30, \) 
90 FORMAT(F20.10) 


WRITE(*,100) 

100 FORMAT(//,2X,’ IS THE LOAD LINEAR OR PARABOLIC (L/P) ? ’,\) 
READ(*,60) RESP 
IF (RESP .BQ. ’L’) THEN 


LD(5) = 1.0 
ELSE 

LD(5) = 2.0 
ENDIF 


ENDIF 


WRITE(*,110) 
110 FORMAT(//,2X,’DO YOU WISH TO MODIFY ANY RING LOADS (Y/N) ? ’,\) 
READ(*,60) RESP 
TF (RESP “EO. (.Y.) EBEN 
WRITE(*,120) 
120 FORMAT(//,2X,’ THE NUMBER OF RING LOADS (INTEGER) = ’\) 
READ(*,130, N 
130 FORMAT(1I3) 
LD(6)=REAL(N) 
K=5 
DO 200 I=1,N 
K=K+2Z 
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WRITE(*,60) ° ° 

WRITE(*,210) ‘MAGNITUDE OF RING LOAD’,I, ’= ’ 
READ(*#,90) LD(KE) 

WRITE(*,210) ‘LOCATION OF RING LOAD ’,I, ’= ’ 
READ(*#,90) LD(K+1) 


200 CONTINUE 
210 FORMAT(2X,A23,13,A2,\) 
ENDIF 
C 


OPEN(7, FILE=’LOAD.DAT’,STATUS=’ NEW’ ) 
DO 219 I=1,50 
WRITE(7,90) LD(T) 
219 CONTINUE 
ENDIF 


WRITE(*#, 220) 
220 FORMAT(//,6X,’ YOU ARE LEAVING THE LOADING MODULE’ ,///) 
END 
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Cc 
Cx 
OF SSS SPS SSS TESST SS SSSTSS SSS SSSSSSSSS SSS SSSSSS SES SSIS SES SSS SS EE SESS SS | 
Cx 
Cx CALL THE SUBROUTINE YLD, WHICH WILL RETURN A LOAD MATRIX 
Cx SUCH THAT YIELDING HAS JUST STARTED TO OCCUR IN THE PLATE 
Cx 
ob SS SS 2S SSS SSS 2S SS 5529535355552 5555225552995335535555 2555595555552 5 5 : 
Cx 
Cc 

CALL YLD(KE, IMAX, LDMAT1) 
Cc 
Cx 
eo SSS 5523955553 555522553555552555525 25525255525 5555253555535222553 8 $ 
Cx 
Cx REPEAT THE PROCESS, AND THEN CALL UPDATE TO REVISE OUR STATE 
Cx MATRIX TO REFLECT THE CONDITION OF THE PLATE AT THE ONSET OF 
C* YIELDING 
Cx 
ee FSS SPSS SSS TSS SSSSSS SSS SSS SSS SSSIS SSS S SESS SIS SESSSESISS SESE SS IESE SSS | 
Cx 

CALL AMATR(BFLAG],BFLAG2,BC,BCPOS, POS, LDMAT1], LHS, AMAT) 

CALL GAUSS(LHS, AMAT, X) 

CALL BCM(POS,BC,X, BCMAT) 

CALL SOLVER(LDMAT1, BCMAT) 

CALL UPDATE(STT1, IMAX,OUT, LDMAT],LL, LDMATO) 

CALL RESUL(M,STT1,COUNT, LDMATO) 
Cc 
Cx 
CE SSTSSSISSS SSS SSS SSSSSSSSESS SSS SSS SSSSSSESSSESSES ESS SESEFESESELSS SLES IS SS 8: 
Ct 
Cx NOW, INPUT AN INCREMENTAL LOAD, AND COMPUTE THE RESULTING 
C* MOMENTS,DEFLECTIONS ETC, USING AS AN INPUT THE PLASTIC MOMENTS 
Cx (IF ANY) AND OTHER DATA FROM THE PREVIOUS LOAD CASE 
Cx 
CE STSSSSSS SSS SESS SSSS SESE SFSFS25955523 952355535295 2359525225595223353522 92s | 
Cx 
Cc 

DO 110 I=1,NN 

LDMAT1(1)=ADJ(5) *LDMATO(1) 
110 CONTINUE 


2) 


116 CONTINUE 


WRITE(*,130) STT(IMAX,9),STT(IMAX, 10) 
130 FORMAT(5X,’MRP IN = °,F20.10,’ MTP IN = °,F20.10) 


aang 


CALL AMATR(BFLAG],BFLAG2,BC,BCPOS, POS, LDMAT1, LHS, AMAT) 
CALL GAUSS(LHS, AMAT, X) 

CALL BCM(POS,BC,X,BCMAT) 

CALL SOLVER(LDMAT1,BCMAT) 
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C 
Cx 


CHAKA AREA RAAT AAR KAKA KER ARAK AAR KERR EKA EAE KAA AKEREEES 


Cx 
Cx 


GO INTO THE SUBROUTINE PLAST TO CALCULATE THE INCREMENTAL PLASTIC 


Cx MOMENTS AT EACH STATION ACROSS THE PLATE WHERE YIELDING HAS 

C* OCCURRED. THE PLASTIC MOMENTS THAT ARE RETURNED BY PLAST ARE 
Cx SQUARED AND ADDED TO THE TOTAL IN THE VARIABLE RMS. RMS IS 

Cx COMPARED TO THE RMS VALUE FROM THE PREVIOUS ITERATION (SAVED AS 


Cx RMS)) | 
Cx ; 
CERKKKAKKKAKAKKKAKAAAKREKATRKKAEKAKEKAAKEKREKAKEKKSEKKKAKEKKEKKEKEKAEKAEKKKK KKK 
Cx 
C 
RMS1=RMS 
RMS=0.0 
IZ=NINT(ADJ(6) ) 
DO 200 I=1,NN 
STTI(1,1)=STT(1,9) 
SrriCl.2)=STT(1,10) 
CALL PLAST(I) 
IF (SET(I,1IZ) .GE. GEOM(6)) THEN 
RMS=RMS+STT(1I,9)**2.0+STT(I,10)*#*2.0 
ENDIF 
200 CONTINUE 
c 
WRITE(*,810) RMS 
810 FORMAT(5X,’RMS = °,F20.10) 
C 
Cx 
Cx THIS WAS THROWN IN TO AVOID A DIVIDE BY ZERO, UNTIL THE PROCESS 
Cx HAS BEEN “PRIMED” 
Cx 
C 
IF (RMS .EG. 0.0) THEN 
RMS=1.0 
ENDIF 
c 
IF (ABS(SQRT(RMS)-SQRT(RMS1))/SQRT(RMS)*100.0 .GT. ADJ(7)) THEN 
DO 815 I=1,NN 
STT(I,9)=(STT(I,9)+STT1(1I,1))/2.0 
STT(1,10)=(STIC1 110) +S72T1 (1,2) )/2.0 
815 CONTINUE 
GOTO 116 
ENDIF 
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Cc 

Cx 

CHEETA ATTA ATES ATA ETAT AE KAKA KARTE KATE AKETAKK EET 
Cx 

om CALL THE SUBROUTINE UPDATE, WHICH WILL UPDATE OUR STATE 

C* MATRIX AND OUR LOADING MATRIX TO REFLECT THE CONDITION OF THE 

Cx PLATE AT THE END OF THE LOADING INCREMENT. TEE PARAMETER 

Cx COUNT KEEPS TRACK OF LOAD INCREMENTS FOR PURPOSES OF DISPLAYING 
Cx RESULTS. WHEN COUNT=10, RESULTS ARE DISPLAYED AND WRITTEN TO 

Cx FILE. (IE EVERY 10 LOAD INCREMENTS ) 


Cx 
SSS SSS SS SS SESS SSS SS SSS SS SSS SS SS SSS SSE S SSS SS SS SSSESSSS SSS SSS PSS SSS 3: 
Cx 
Cc 
CALL UPDATE(STT1,IMAX,OUT,LDMAT1,LL,LDMATO) 
CALL RESUL(M,STT1,COUNT, LDMATO) 
IF (COUNT .EQ. 10.0) THEN 
COUNT=0.0 
ENDIF 
COUNT=COUNT+1.0 
M=M+] 
WRITE(*,117) M 
117 FORMAT(5X,’ LOAD INCR = ’,1I3) 
& 
es 


CEEEKKEEEEEKKEEEEKEK ETEK ERE EKKK STAKE KEKEKAKK ERE RK KEKE KKK KKK 
Cx 
Cx CHECK TO SEE IF EXCEEDING ULTIMATE STRAIN. IF SO, KICK 
Cx OUT OF THE PROCEDURE. IF NOT, WE SHOULD INCREMENT ON UP AGAIN 
Cx BY USING AS OUR INPUT THE INCREMENTAL LOAD MATRIX. 
Cx 
CEEEKKEKEKEEKEKEKEEEKAAEKK ATER EK EATER EK AEK KATE KAKKAKS 
Cx 

IF (OUT .NE. 1) THEN 


GOTO 116 

ENDIF 
C 

CALL RESUL(M,STT1,COUNT, LDMATO) 
C 

CLOSE(15) 
C 

END 
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Cc 
Cx 
CEEEEEESEEEE EEE TEE KEKTEKEE TEAK KEKE EERE RAKE KEKE EEE EERE KEKE KE 
cx 


Cx SUBROUTINE INTEG 
Cx 
Cx THIS SUBROUTINE IS CALLED BY THE PROGRAM MAIN TO CREATE THE 
Cx ARRAY LDMAT. THIS ARRAY CONTAINS THE VALUES OF THE TOTAL LOAD 
Cx FOR EACH SEGMENT ALONG THE PLATE RADIUS. IT IS CONSTRUCTED 
cx BY USING THE AVERAGE OF THE TWO ENDPOINTS FOR DISTRIBUTED 
Cx LOADS AND ADDING TO THIS VALUE THE MAGNITUDE OF ANY RING LOADS 
Cx THAT ARE SITUATED WITHIN THE SEGMENT. 
Cx 
CEEXKXKKKKEKESESEEEEKEKEKEKEEE EE EKEKEEKEKAE EEE EERE KEE EKER EKEKEKEKSE 
Cx 
C 
SUBROUTINE INTEG(LD,LDMAT) 
C 
IMPLICIT REAL*8 (A-H,0-Z) 
REAL*8 LDMAT(50),LD(50),NU 
INTEGER N 
SINCLUDE: ’COMMON. FOR’ 
C 
A=GEOM(1) 
B=GEOM(2) 
NU=GEOM(3) 
D=GEOM(4) 
C 
K=NINT(ADJ(2)) 
DEL=(A-B)/ADJ(2) 
R=B 
N=NINT(6.0+2.0¥8LD(6) ) 
Cc 
Cr 


CEE KEKKEEKEKEKEEE KE KRER ERA EAE KEKE KEES 
cx 


Cx CHECK TO SEE IF THERE IS A DISTRIBUTED LOAD. IF NOT, PROCEED 
C* TO THE SECTION FOR RING LOADS 
Cx 


CEEKKKKXEEKEKEKEEEEKAEE EEE KEKE RAKE EKER CER EEE KEKE ERE EEE EEE EEE 
om 
Cc 
IF (LD(1) .NE. 0.0 .OR. LD(3) .NE. 0.0) THEN 

C 

R=B 

SLP=(LD(1)-LD(3))/(LD(2)-LD(4)) 

DO 100 I=1,K 


0 





C 
Cx 
1 2S 2S SS SSS 2 SSS SLSSSSP PSPS LL SS SLL SLPS SPSS S LES LSSE LES LP SSPP LSE LE SESS SSS 2 & @ 


Cx SECTION FOR LINEARLY DISTRIBUTED LOADS 


S222 FS SSS SSP SPSL? LS LP LSP LSP LESPS LES SSPP LS SS SE LLLPSSSLSLLSSSLSL SSE LES SSDS fF © : 
Cx 


C 
IF (LD(5) .EQ. 1.0 ) THEN 
C 
R=R+DEL 
C 
IF (R-LD(2) .GT. DEL/10.0) THEN 
LDMAT(1)=0.0 
C 
ELSE IF (R .GT. LD(4)) THEN 
LDMAT(1I)=( LP*(R-LD(4)-DEL/2.0)+LD(3))*DEL 
C 
ELSE 
C 
LOMAT(1I)=0.0 
C 
ENDIF 
C 
Cx 


OF SSP SP SSSSSSSSSSSSSSSSSS LSP LIL SS SLSPSPLLLIS LESS SLL SPSS LLP SSS 2 SSS Sf SF ® F 
Cx 

Cx SECTION FOR PARABOLIC LOADS (Y=X*x2) 

Cx 

Of SSS SSS SSS SEP SPPSSPSSSCSS SSP L LS SSP LIES LL SSS LPLS LP SELES S PPS SL LS SPSS SSS ® F | 
Cx 


C 
ELSE IF (LD(5) .EQ. 2.0 ) THEN 
C 
R=R+DEL 
P=(LD(2)-LD(4) )**2.0/4.0/(LD(1)-LD(3)) 
IF (LD(1) .GT. LD(3)) THEN 
Y=(R-DEL/2.0-LD(4) ) **2.0/4.0/P+LD(3) 
ELSE 
Y=-(LD(2)-R+DEL/2.0)**2.0/4.0/P+LD(4) 
ENDIF 
C 
IF (R-LD(2) .GT. DEL/10.0) THEN 
LDMAT(1I)=0.0 
C 
ELSE IF (R .GT. LD(4)) THEN 
LDMAT(I)=Y*DEL 
C 
ELSE 
C 
LDMAT(1I)=0.0 
C 
ENDIF 
ENDIF 


ited 





¢ 

100 CONTINUE 

ENDIF 

c 
Cx 
CEEXKEEEKEEKKEKKEEKRAK KAKA KAKA EEE AEE AEE ERE EEK EEE KEE EEE EES 
cx 
Cx SECTION FOR ADDING IN RING LOADS 
Cx 
Cex EKEKEKEK EKA EKA EEE EAE KEKE KEKE KEKE 
Cx 


Cc 
J=5 
IF (N .NE. 0) THEN 
DO 300 I=7,6+NINT(LD(6) ) 
J=J+2 
L=INT((LD(J+1)-B)/DEL)+1 
LDMAT(L)=LDMAT(L)+LD(J) 
300 CONTINUE 
ENDIF 
C 
c 
RETURN 
END 
c 
Cx 


CEOEEKEKEKEAEKEKEKKEKEKR EEK EK KAKEEEEKEKE KE ERE KEKEKEKEEEEEKEKEEEEKEKSE 
Cx 

Cx SUBROUTINE AMATR 

Cx 

Cx THIS SUBROUTINE IS CALLED BY THE PROGRAM MAIN TO CREATE THE 

Cx 4X4 A MATRIX, WHICH IS SOLVED USING GAUSS ELIMINATION TO 

C* DETERMINE THE UNKNOWN BOUNDARY CONDITIONS. 


Cx 
CeKEKEKERKEKEKEKEKEKARKEEEKEEEKEKEKEK EAE EEE EEE KEE EK EKER EE EEKEKEKEKSE 
Cx 
Cc 
SUBROUTINE AMATR(BFLAG1,BFLAG2,BC,BCPOS, POS, LDMAT1, LHS, AMAT) 
c 
IMPLICIT REAL*¥8 (A-H,0-Z) 
CHARACTER*1 BFLAG1(4), BFLAG2(4) 
REAL*8 BC(4),LHS(4),INC,NU,GRNMT1(4,6),GRNMT2(4,6), 
+ AMAT(4,4),LDMAT(4),LDTOT(4), LDMAT1 (50) ,MPRTOT(4) ,MPTTOT(4), 
- STTAV(5S07,2) 
INTEGER POS(8),BCPOS(8) ,M,N 
C 
$INCLUDE: ’COMMON. FOR’ 
c 


A=GEOM(1]) 


eZ 





B=GEOM(2) 


NU=GEOM(3) 
D=GEOM(4) 
INC=ADJ(1) 
C 
DO 101 I=1,4 
LHS(1I)=0.0 
IF (I .EQ. 1 .OR. I .BQ. 3) THEN 
RO=A-INC 
ELSE 
RO=B+INC 
ENDIF 
C 
R=A 
CALL GRNFNC(R,RO,NU,D,GRNMT1) 
R=B 
CALL GRNFNC(R,RO,NU,D,GRNMT2) 
% 
Cx 
CEEEEEEEEKEEKEEEEEKEKREKEEEEEEE EKER AEA EEE AAA EKEAKEKA KAKA EKELES 
Cx 
Cx REDUCE THE TERMS ORIGINALLY ON THE RHS TO A 4X4 MATRIX 


Cx COEFFICIENTS (AMAT) AND A 4X1 MATRIX OF NON-COEFFICIENT 

Cx TERMS (LHS) 

Cx 
CEEEXKKEKEKKEKEEAEA EEK EAA KARATE ERE RATA EATER AEA EKER ATES 
Cx 

Cc 


thy (Tl e021 OR. 1 wEOQs 2) THEN 


IF (BFLAG1(J) .NE. *U’) THEN 
K=K+] 
LHS(I)=GRNMT1(J,1)*BC(K)*A*(-1.0)4%%(J-1)+LHS(1I) 
ELSE 
L=L+} 
AMAT(I,L)=GRNMT1(J,1)*A%(-1.0)**%(J) 
ENDIF 
73 FORMAT(2X,F10.5) 
90 CONTINUE 


DO 100 J=1,4 
IF (BFLAG2(J) .NE. ’°U’) THEN 
K=K+] 
LHS (I)=GRNMT2(J,1)*BC(K) *B*(-1) *4*J+LHS(1) 
ELSE 
L=L+1 
AMAT(1,L)=GRNMT2(J,1)*B*x(-1)**(J-1) 
ENDIF 
100 CONTINUE 


ELSE 
K=0 


es 





L=0 
DO 110 J=1,4 
IF (BFLAGI(J) .NE. ’U’) THEN 
K=K+1 
LHS(1)=GRNMT1(J,2)*BC(K) #A¥(-1.0)*%*(J-1)+LHS(1) 
ELSE 
L=L+1 
AMAT(1I,L)=GRNMT1(J,2) *A*(-1.0) *8(J) 
ENDIF 
110 CONTINUE 


DO 115 J=1,4 
IF (BFLAG2(J) .NE. °U’) THEN 
K=K+] 
LHS(I)=GRNMT2(J,2)*BC(K)*B2*(-1)**(J)+LHS(I) 
ELSE 
L=L+1 
AMAT(I,L)=GRNMT2(J,2)*B*(-1)**«(J-1) 
ENDIF 
115 CONTINUE 
ENDIF 
Cc 
Cx 
CEKEKEEKEKESEKEKEEKAEKEKREKE KEKE EEE KEK ERE KER EERE KEKE EKA EEE EEE ERE EEK 
Cx 
Cx MOVE THE TERMS ORIGINALLY ON THE LHS TO THE RHS, FORM A NEW AMAT 
C* MATRIX AND A NEW LHS MATRIX WHICH CORRESPOND TO THE A AND B MATRICES 
Cx OF THE SYSTEM AX=B 
Cx 
CExKEEKKKEEEKEEKEKEKEAREKAAREKEKEKR KEKE KAKA EEA ERK SET ARERR EKER EKEKEKEKKE 
Cx 
IF (BFLAG1(4) .NE. °U’ . AND. I .EQ. 1) THEN 
LHS(1)=LHS(1)+RO*BC(BCPOS(4) ) 
ELSE IF (BFLAG1(4) .BQ. ’°U’ .AND. I .EQ. 1) THEN 
AMAT(I,POS(4))=AMAT(I,POS(4))-RO 


ENDIF 
C 
IF (BFLAG2(4) .NE. ’U’ .AND. I .EQ. 2) THEN 
LHS(2)=LHS(2)+RO*BC(BCPOS(8) ) 
ELSE IF (BFLAG2(4) .EQ. ’U’ .AND. I .EQ. 2) THEN 
AMAT(1I,POS(8))=AMAT(I, POS(8))-RO 
ENDIF 
Cc 


IF (BFLAGI(3) .NE. ’°U’ .AND. I .EQ. 3) THEN 
LHS(3)=LHS(3)-RO*BC(BCPOS(3)) 

ENDIF 

IF (BFLAG1(4) .NE. °U’ .AND. JI .EQ. 3) THEN 
LHS(3)=LHS(3)+BC(BCPOS( 4) ) 

ENDIF 

IF (BFLAGI(3) .EQ@. ’*°U’ .AND. I .BQ. 3) THEN 
AMAT(I,POS(3))=AMAT(I,POS(3))+RO 

ENDIF 

IF (BFLAG1(4) .EQ. °U’ .AND. I .EQ. 3) THEN 
AMAT(I,POS(4)}=AMAT(I,POS(4)}-1.0 

ENDIF 
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IF (BFLAG2(3) .NE. 'U’ .AND. I .EQ. 4) THEN 
LHS (4)=LHS(4)-RO*BC(BCPOS(7) ) 

ENDIF 

IF (BFLAG2(4) .NE. 'U’ .AND. I .EQ. 4) THEN 
LHS (4) =LHS(4)+BC(BCPOS(8) ) 

ENDIF 

The CBELAGZ( 3) .EQ. “U° .AND. I .EQ. 4) THEN 
AMAT(I,POS(7))=AMAT(I,POS(7))+RO 

ENDIF 

IF (BFLAG2(4) .EQ. °U’ .AND. I .EQ. 4) THEN 
AMAT(I,POS(8))=AMAT(I,POS(8))-1.0 


ENDIF 
Cc 
101 CONTINUE 
c 
Cx 
SSS SCS SS SSS SESS SSS SSS SS SS SSS SCS SSCS SS SESSS CSS SSS SCS SS SSCS SSS SESS: 
Cr 
Cx FINALLY, ADJUST THE LHS TO ACCOUNT FOR THE EXTERNAL LOAD 
C* AND PLASTICITY BY CALCULATING AND SUBTRACTING APPROPRIATE 
C* TERMS. 
Cx 
SSS SSS SSCS SSCS SSS SSSSCS SSCS SSS SSCS SS SSS SSS SES SSCS S SSCS SSE SSS SSS SESS: 
Cx 
C 
K=NINT(ADJ(2) ) 
DEL=(A-B)/ADJ(2) 
C 
DO 500 I=1,Kk F 

R=B-DEL/2.0+REAL( I) *DEL 

RO=A-ADJ(1) 

STIAVCL, F)=(STI(1, S)4STIC1+1,9))72.0 

SITAV Cr. 2)=CSTTCt, 10) 7stT(i+1,10))/2.0 

CALL GRNFNC(R,RO,NU,D,GRNMT1) 
LDTOT(1)=LDTOT(1)+LDMAT1(1I) *GRNMT1(1,1)#R 
MPRTOT(1)=MPRTOT(1)+STTAV(I,1)*(-GRNMT1(3,1)/D+ 

+ NU/R*GRNMT1(2,1))*DEL*R 
MPTTOT(1)=MPTTOT(1)+STTAV(I,2)*(-GRNMT1(2,1)*DEL) 
LDTOT(3)=LDTOT(3)+LDMAT1(1) *GRNMT1(1,2)%*R 
MPRTOT(3)=MPRTOT(3)+STTAV(I,1)*(-GRNMT1(3,2)/D+ 

+ NU/R*¥GRNMT1(2,2))*DEL*R 
MPTTOT(3)=MPTTOT(3)+STTAV(I,2)*(-GRNMT1(2,2)*DEL) 

RO=B+ADJ(1) 

CALL GRNFNC(R,RO,NU,D,GRNMT1) 
LDTOT(2)=LDTOT(2)+LDMAT1(1)*GRNMT1(1,1)4R 
MPRTOT(2)=MPRTOT(2)+STTAV(1I,1)*(-GRNMT1(3,1)/D+ 

+ NU/R*GRNMT1(2,1))*DEL*R 
MPTTOT(2)=MPTTOT(2)+STTAV(I,2)*(-GRNMT1(2,1)*DEL) 
LDTOT(4)=LDTOT(4)+LDMAT1(1)*GRNMT1(1,2)#R 
MPRTOT(4)=MPRTOT(4)+STTAV(1I,1)*(-GRNMT1(3,2)/D+ 

+ NU/R*GRNMT1(2,2))*DEL*R 
MPTTOT(4)=MPTTOT(4)+STTAV(I,2)*(-GRNMT1(2,2)xDEL) 

500 CONTINUE 
C 
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DO 105 1=1,4 
LHS (I)=LHS(I)-LDTOT(I)+MPRTOT(I)+MPTTOT(I) 
105 CONTINUE 


& 
DO 1000 I=1,4 
LDTOT(I)=0.0 
MPRTOT(I)=0.0 
MPTTOT(I)=0.0 
1000 CONTINUE 
c 
RETURN 
END 
Cc 
Cx 
Ob SSS SSS SSS SSS SSS SSS SSS S STS S SSS S SSS SESS SS SSS SSS SSS S SSS SS SESS StS s; 
Cx 
cx SUBROUTINE GRNFNC 
cr 
Cx THIS SUBROUTINE CALCULATES THE VALUES OF DEFLECTION, SLOPE, 


Cx MOMENT AND SHEAR AND ASSOCIATED DERIVATIVES FOR THE GREEN’S 
Cx FUNCTION, AND RETURNS THESE VALUES TO THE MAIN PROGRAM. TO 
Cx INDICATE DERIVATIVES OF PARAMETERS, THE NUMBER 1, 2, OR 3 IS 
Cx TACKED ON TO THE END OF THE VARIABLE NAME. HENCE, THE FIRST 
C* DERIVATIVE OF LS WRT RO IS L392 
Cr 
Cx 
SS SSS SSS SESS SSS SSCS SS SSS SSCS SESS SESE SESE SESS SESE SOS S SSS e Ses eeseess: 
cx 
C 

SUBROUTINE GRNFNC(R,RO,NU,D,GRNMAT) 
© 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 NU,GRNMAT(4,6),1L3,L31,L32,L33,L9,L91,L92,1L93 
S$INCLUDE: ’COMMON. FOR’ 
C 

w=1.0 

A=GEOM(1)+ADJ(4) 

B=GEOM(2)-ADJ(4) 

THB=0.0 

MB=0.0 

QB=0.0 


IF (R .GT. RO) THEN 
RRO=1.0 

ELSE 
RRO=0.0 

ENDIF 
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C 

Ct 

CEEKETEKEKETA TATE TAAAAAAKAATE KATA KKAAAKTKT TAKER EKER AA RAAT 
Cx 


cx CALCULATE PLATE CONSTANTS AND PLATE FUNCTIONS, AND THEIR 
Cx DERIVATIVES 
Ct 


CEEEKEREEEK AEE AERA AAA AEST AE AEKATA KAA AAAKE RAAT AKA AAA AKER KE 

Cx 

Cc 5 
L3=RO/4.0/A*(((RO/A) #2. 0+1.0) *DLOG(A/RO)+(RO/A) *%2.0-1.0) 
L31=1.0/4.0/A*(DLOG(A/RO) &(3.08R0*%2.0/A*42.0+1.0)+2.0% 

+ (BO**2.0/A%%2.0-1.0)) 
L32=1.0/4.0/A#(6.0RO/A*%2.0*#DLOG(A/RO)+RO/A#%2.0-1.0/R0) 
L33=1.0/4.0/A%(6.0/A##2.0%(DLOG(A/RO)-1.0)+1.0/A%42.0+1.0/RO 

+ %*%*2.0) 


LO=RO/A¥((1.0+NU)/2.02DLOG(A/RO)+(1.0-NU) /4.0%(1.0-(RO/A) **2.0)) 
L91=1.0/4.0/A*(2.04(1.0+NU) #DLOG(A/RO)-(3.0#NU+1.0)-3.0%(1.0-NU) 
+ *RO#%2.0/A*4#2.0) 
L92=-(1.0+NU)/2.0/A4/RO-3.0#RO#(1.0-NU)/2.0/A*¥¥3.0 
L93=(1.0+NU)/2.0/A/RO##2.0-3.0#(1.0-NU)/2.0/A*¥3.0 


Cl=(1.0+NU)/2.0*#B/A*DLOG(A/B)+(1.0-NU)/4.0%(A/B-B/A) 
C7=1.0/2.0%(1.0-NU*%2.0) &(A/B-B/A) 


G3=RO/4.0/R%(((RO/R) #*2.0+1.0) SD LOG(R/RO)+(RO/R) **2.0-1.0) RRO 
G31=1.0/4.0/R%*(3.0%RO*#2.0/R422.02DL0G(R/RO)+2.04RO042%2.0/R422.0+ 

+ DLOG(R/RO)-2.0) *RRO 
G32=1.0/4.0/R%*(6.0*RO/R*22.0*DL0G(R/RO)+RO/R*¥*2.0-1.0/RO) *RRO 
G33=1.0/4.0/R%(6.0/Rx*2.0%(DLOG(R/RO)-1.0)+1.0/R¥22.0+1.0/RO%%2.0) 
+  xXRRO 


G6=RO/4.0/R*((RO/R) *22.0-1.0+2.0DLOG(R/RO) ) RRO 
G61=1.0/4.0/R*(3.0%(ROs%2.0/R¥42.0-1.0)+2.0DLOG(R/RO) ) RRO 
G62=(3.0*%R0O/2.0/R2%*3.0-1.0/2.0/RO/R) xRRO 
G63=(3.0/2.0/BR*x*x3.0+1.0/2.0/R/RO*#%x2.0) *RRO 


G9=RO/R*((1.0+NU)/2.0*DLOG(R/RO)+(1.0-NU) /4.04(1.0-(RO/R) *%2.0) ) 
+ *RRO 

G91=1.0/4.0/R*«(2.0%(1.0+NU) x(DLOG(R/RO)-1.0)+(1.0-NU)-3.0% 

+ (1.0-NU) #RO*#%2.0/R422.0) RRO 
G92=-1.0/4.0/R*(2.04#(1.0+NU) /RO+6. OFRO/R*¥4*2.04(1.0-NU) ) *RRO 
G93=1.0/4.0/R«(2.0%(1.0+NU) /RO**#2.0-6.0/Rx%2.0%(1.0-NU) ) *RRO 


Fl=(1.0+NU)/2.0%B/R&DLOG(R/B)+(1.0-NU)/4.0%(R/B-B/R) 
F2=1.0/4.0%(1.0-(B/R) **2.0%(1.0+2.0*8DLOG(R/B) ) ) 
F3=B/4.0/R*(((B/R)**2.0+1.0) *DLOG(R/B)+(B/R)**2.0-1.0) 
F4=1.0/2.0%((1.0+NU) *B/R+(1.0-NU) ¥R/B) 
F5=1.0/2.0%(1.0-(B/R) **2.0) 

F6=B/4.0/R*((B/R) **2.0-1.0+2.08DLOG(R/B) ) 
F7=1.0/2.0%(1.0-NU%%&2.0)%*(R/B-B/R) 
F8=1.0/2.0%(1.0+NU+(1.0-NU)*(B/R) **2.0) 
FQ9=B/R*¥((1.0+NU)/2.0*%DLOG(R/B)+(1.0-NU)/4.0%(1.0-(B/R)**2.0) ) 


La 








C 
Ct 


Oe FSP SSPE LISPIPLISILS SPSS SPS ILS SPSLISPE SSPE SPPSESSESIS SS SSSSSSSSSSSSSSS SS SSS | 


Cx 
Ct 
Ct 
Cx 
Cx 
Cx 
Cx 
Cx 
Cx 
Cx 
Cx 
Cx 
Cx 
Cx 
Ct 
C* 
Cx 
Ct 
Ct 
Ct 
C* 
Ct 


CALCULATE THE DEFLECTION, SLOPE, MOMENT AND SHEAR, AND 
ASSOCIATED DERIVATIVES, AND STORE IN THE MATRIX GRNMAT. 
THIS MATRIX IS A 4 BY 4 MATRIX, WHERE ROW POSITIONS ARE AS 
FOLLOWS 


1 DEFLECTION 
2 SLOPE 

3 MOMENT 

4 SHEAR 


AND COLUMN POSITIONS ARE AS FOLLOWS 


1 GREEN’S FUNCTION (NO DERIVATIVES) 

2 1ST DERIVATIVE OF GREEN’S FUNCTION WRT RO 
3 2ND DERIVATIVE OF GREEN’S FUNCTION WRT RO 
4 3RD DERIVATIVE OF GREEN’S FUNCTION WRT RO 


FOR EXAMPLE, THE VALUE RETURNED BY GRNMAT(2,3) REPRESENTS 
THE SECOND DERIVATIVE OF THE SLOPE FOR THE GREEN’S FUNCTION 
WRT RO 


2 FSS PPPS SPL LIPLISSLLSLSLISSSSSE PISS LPS LISLELSPESLPSS PES LLP SPSL PLL SSS SS SS ® ® ® 


Ct 
C 


YB=-WtA*%3.0/D*(C1*L9/C7-L3) 
THB=W*A%42.0/D/C7*#LI 
GRNMAT(1,1)=-(YB+THB*R*F1-WeR**3.0/D%G3) 
GRNMAT(2,1)=THB*F4-W*R*%2.0/D4G6 
GRNMAT(3,1)=THB*D/R*F7-WtR*G9 
GRNMAT(4,1)=W*RO/R*RRO 


YBl=-W*A¥4%3.0/D4(C1*L91/C7-L31) 
THBl=W*A¥4%2.0/D/C74#L91 
GRNMAT(1,2)=-(YBI1+THBI*R*F1-WtR4¥43.0/D4G31) 
GRNMAT(2,2)=THB1*F4-W*R4*4*2.0/D*G61 
GRNMAT(3,2)=THB1*D/R*F7-W*R*G91] 
GRNMAT(4,2)=W/R*BRO 


YB2=-W*tA¥43.0/D%(C1/C7*L92-L32) 
THB2=WHAX*2.0/D/C7*LI92 

GRNMAT(1,3)=-( YB2+R#F1*THB2-W#R*¥43.0/D*G32) 
GRNMAT(2,3)=F4*THB2-W*xR*42.0/D*G62 
GRNMAT(3,3)=D#F7/R*THB2-W*R*¥G92 
GRNMAT(4,3)=0 


YB3=-W*A*%3.0/D%(C1/C7*L93-L33) 
THB3=W*A¥%2.0/D/C7#L93 
GRNMAT(1,4)=-(YB3+R*F1*THB3-W*R*¥*3.0/D%G33; 
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GRNMAT(2,4)=F4xTHB3-W*eh¥42.0/D%G63 
GRNMAT(3,4)=D*eF7/R*THB3-W*R*G93 
GRNMAT(4,4)=0 


C 
RETUEN 
END 
C 
cx 
CEEKEEKEEEEKEKEKREKEKEKEKEEEEKEKKEKEKEREKEKEKREKKEREEEERKEKEKEEEKREEEEREKREKEE 
Cx 
Cx SUBROUTINE GAUSS 
Cx 
Cx THIS SUBROUTINE USES GAUSS ELIMINATION TO DETERMINE THE 


Cx UNKNOWN BOUNDARY CONDITIONS. THE MATRICES A AND B OF THE 
Cx EQUATION AX=B ARE PROVIDED BY THE SUBROUTINE AMATR. THIS 
Cx SUBROUTINE THEN RETURNS THE MATRIX X OF UNKNOWN BOUNDARY 
Cx CONDITIONS. THIS SUBROUTINE WAS TAKEN FROM THE BOOK 
Ct “NUMERICAL ANALYSIS” BY L.W. JOHNSON AND R.D. RIESS. 


Cx 
CEEEEEEKEKKEKEKEEEREKEKKEEKKEKREKEREKEKERERKEKEKKKEREKAERERKEEKEEKEEEREREREKEKEKESE 
Cx 
c 
SUBROUTINE GAUSS(B,A,X) 
C 
IMPLICIT REAL*8 (A-H,0-Z) 
CHARACTER*22 BCLBL(8) 
REAL*8 A(4,4),B(4),X(4),BCMAT(8),BC(4) 
DIMENSION AUG(50,51) 
INTEGER POS(8),BCPOS(8),M,N 
Cc 
NM1=3 
NP=5 
N=4 
C 
C SET UP THE AUGMENTED MATRIX FOR AX=B 
C 
DO 2 I=1,N 
DO 1 J=1,N 
AUG(I,J)=A(I,J) 
if CONTINUE 


AUG(I,NP)=B(I) 
2 CONTINUE 


THE OUTER LOOP USES ELEMENTARY ROW OPERATIONS TO TRANSFORM 
THE AUGMENTED MATRIX TO ECHELON FORM 


YPQaAD 
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oO 


DO 8 I=1,NM1 


SEARCH FOR THE LARGEST ENTRY IN COLUMN I, ROWS I THROUGH N 
IPIVOT IS THE ROW INDEX OF THE LARGEST ENTRY 


AAA 


PIVOT=0.0 : 
DOs 30 =a, N 
TEMP=ABS(AUG(J,1)) 
IF(PIVOT.GE.TEMP) GO TO 3 
PIVOT=TEMP 
IPIVOT=J 
3 CONTINUE 
IF(PIVOT.EQ.0.0) GO TO 13 
IF(IPIVOT.EQ.I) GO TO 5 
C 
C INTERCHANGE ROW I AND ROW IPIVOT 
Cc 
DO 4 K=I,NP 
TEMP=AUG(I,K) 
AUG(I,K)=AUG(IPIVOT,&) 
AUG(IPIVOT,K)=TEMP 
4 CONTINUE 
Cc 
C ZERO ENTRIES (I+1),(i+2),...,(N,1) IN THE AUGMENTED MATRIX 
Cc 
5 TPl=1+1 
DO 7 K=IP1,N 
Q=-AUG(K,1I)/AUG(I,I) 
AUG(K, 1)=0.0 
DO 6 J=IP1,NP 
AUG(K,J)=Q*AUG(I,J)+AUG(K, J) 
- CONTINUE 
i CONTINUE 
8 CONTINUE 
IF(AUG(N,N).EQ.0.0) GO TO 13 
c 
C BACKSOLVE TO OBTAIN A SOLUTION TO AX=B 
c 
X(N)=AUG(N,NP)/AUG(N,N) 
DO 10 K=1,NMl 
e=0 0 
DO 9 J=1,K 
Q@=Q+AUG(N-K,NP-J) #®X(NP-J) 
9 CONTINUE 
X (N-K)=(AUG(N-K, NP)-Q) /AUG(N-K,N-K) 
10 CONTINUE 


CALCULATE THE NORM OF THE RESIDUAL VECTOR, B-AX. 
SET IERROR=1 AND RETURN 


aaAANgn 


RSQ=0.0 
DO 12 I=1,N 
Q=0.0 
DO 11 J=1,N 
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Q=Q+A(I,J)¥X(J) 
io} CONTINUE 
REST=B(I)-@Q 
RMAG=ABS(RBSI) 
RSQ=RSQ+RMAG*¥2 
2 CONTINUE 
RNORM=SQRT(RSQ) 
IERROR=1] 


ABNORMAL RETURN --- REDUCTION TO ECHELON FORM PRODUCES A ZERO 
ENTRY ON THE DIAGONAL. THE MATRIX A MAY BE SINGULAR. 


13 IERROR=2 


QO OaAa an 


RETURN 
END 


C 

Cx 

SS SSSSSSSLSPSSSSSSSSSSSESSSSSSSSSSSSSSSSSSSPSSSS SSS SSS SSSSSSSSSS Ss Ss ese: 
Cx 

Cx SUBROUTINE BCM 

Cx 

Cx THIS SUBROUTINE TAKES THE KNOWN BOUNDARY CONDITIONS STORED 

Cx BY THE PROGRAM INPUT IN THE FILE BCM.DAT AS WELL AS THE 

Cx "UNKNOWN" BOUNDARY CONDITIONS RETURNED BY THE SUBROUTINE GAUSS 
C* AND CREATES THE 8X1 MATRIX OF BOUNDARY CONDITIONS BCMAT. 

Cx 

CEEKEKKKEKKAEEKEKE EKER KEKE EAE 
Cx 


C 
SUBROUTINE BCM(POS,BC,X,BCMAT) 
C 
IMPLICIT REAL*¥8 (A-H,0-Z)C 
REAL*8 BC(4),X(4),BCMAT(8) 
INTEGER POS(8) 
c 
K=0 
L=0 
DO 120 I1=1,8 
IF (POS(I) .GT. 0) THEN 
K=K+] 
BCMAT(1)=X(K) 
ELSE 
Laie) 
BCMAT(I)=BC(L) 
ENDIF 
120 CONTINUE 
RETURN 


END 
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C 
Ce 
2 FSS SSS SSSSSSSSESSS ES SS SSS SSSSS SPSS SSSSSESSS SSS SSS SSS SSS SS ES SPSS SF FF | 


Cx SUBROUTINE SOLVER 


Cx THIS SUBROUTINE SOLVES TAKES THE COMPLETE SET OF BOUNDARY 

Cx CONDITIONS SUPPLIED BY THE SUBROUTINE BCM AND SOLVES THE PLATE 
Cx BENDING PROBLEM ACROSS THE PLATE WIDTH. IMPORTANT PARAMETERS 

Ct SUCH AS DEFLECTION, SLOPE ETC ARE STORED IN THE ARRAY STT. 

Cx 

OSS SSESESSSPSSSSSSSSSSSSIESSSSSLLSSPSESSLSESSPESESS SESS SSESESS SILLS SSE FS ® ® & F | 


Cx 
C 
SUBROUTINE SOLVER(XMT, BCMAT) 
Cc 
IMPLICIT REAL*¥8 (A-H,0-Z) 
C 


REAL*8 GRNMT1(4,6),GRNMT2(4,6),BCMAT(8),DT(50,4),SM(50,4),NU, 
+ XMT(50),STTAV(S5O, 2) 
INTEGER M,N 

$ INCLUDE: COMMON. FOR’ 
C 
Cx 
OP PSS SSS SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSPSFSSLSSSSSES SSS SS SSS SSS SSDS SSF | 
or 
Cx REZERO THE MATRICES DT AND ST 
cx 
OF SS SPP SSSSSSSSSSSSSSSSSSSSSS SPSS SS SSSSSSSSS SSSI PSSESSSSSSSS SS SSS SSDS SS | 
Cx 


C 
DO 10 I=1,50 
DO 5 J=1,4 
SM(I,J)=0.0 
DT(1I,J)=0.0 
S CONTINUE 
10 CONTINUE 
C 
A=GEOM( 1) 
B=GEOM( 2) 
NU=GEOM(3) 
D=GEOM(4) 
C 
C 
Cc 
RO=B 


DEL=(A-B)/ADJ(2) 
N=INT(ADJ(2)) 


Lge 








C 
Cx 


CHERKEAATETAARAARARTAAA TATRA KKARA RAKE RETA RETAKE REET AEEKKE 


Cx 
Cx 
Cx 
Cx 
Cx 
Cx 
Cx 
Cx 
Cxx* 
Cx 
Cc 


C 


140 
150 


LS5 
C 
& 
Cr 
Cx*xx 
Cx 
Cx 
Cx 
Cr 
Cr 
Cx 
Cwxx 
Cx 
C 


C 


THIS LOOP CALCULATES THE LOAD TERM AND THE PLASTIC MOMENT 
TERM FOR EACH EQUATION. I REFERS TO THE RADIAL STATION AND J 
REFERS TO THE EQUATION, WHERE THE FIRST EQUATION IS FOR 
DEFLECTION, THE SECOND FOR DW/DR, THE THIRD FOR D”“2W/DR“2 AND 
THE FOURTH FOR D“3W/DR*3. THE DERIVATIVES ARE THEN USED TO 
CALCULATE SLOPES, MOMENTS AND SHEARS 


RERKTAAEERETLEKAAE TRAE TAE TATE ETAKA EAA RETAKE REACT AACE ERATE 


DO 155 I=1,N+l 


IF (I .£Q. 1) THEN 
RO=B+ADJ(1) 
ENDIF 
TEOGl - EO. (N+i)) THEN 
RO=A-ADJ(1) 
ENDIF 
DO 150 M=1,NINT(ADJ(2)) 
R=B-DEL/2.0+REAL(M) sDEL 
CALL GRNFNC(R,RO,NU,D,GRNMT1) 
STTAV(M,1)=(STT(M,9)+STT(M+1,9))/2.0 
STTAV(M,2)=(STT(M,10)+STT(M+1,10))/2.0 
DO 140 J=1,4 
DT(I,J)=DT(I,J)+(XMT(M) *GRNMT1(1,J)*R-DEL&*STTAV(M,1)*% 
+ (-GRNMT1(3,J)/D+NU/R&*GRNMT1(2,J))4R 
* +DEL&STTAV(M, 2) *GRNMT1(2,J)) 
CONTINUE 
CONTINUE 
Tes(T EQ. 1) THEN 
RO=BtDEL 
ELS§ 
RO=RO+DEL 
ENDIF 
CONTINUE 


REKKKKA AAA AAR KAKA AA RAAAAA AAA RATER KAA ARA EKA KAKA ERAT AAR TE 
AT EACH INCREMENT OF WIDTH DEL ALONG THE PLATE RADIUS, CALCULATE 

THE DEFLECTION, SLOPE, MOMENT AND SHEAR, WHICH ARE STORED IN 

STT (EACH COLUMN REPRESENTS RO, Y, TH, M, @, Y’’, Y’’’ CORRES- 

PONDING TO EACH STATION I ACROSS THE RADIUS} 


XERKEEEATETAEAEEEA EET EKER KEE EAE ATE EAE ETAT EEK TATE IE 


DO 200 I=1,N+l 


Pea] 2. 2G > 1) THEN 


eZ 








+ 


170 


180 
C 
Cx 


RO=B+ADJ(1) 

ENDIF 

IF (I .BQ. (N+1l1)) THEN 
RO=A-ADJ(1) 

ENDIF 


R=A 

CALL GRNFNC(R,RO,NU,D,GRNMT1) 

R=B 

CALL GRNFNC(R,RO,NU,D,GRNMT2) 

DO 180 J=1,4 
DO 170 K=1,4 

SM(I,J)=A*GRNMT1(K,J)*BCMAT(K)*(-1.0) **K- 
B¥GRNMT2(K,J)*BCMAT(K+4)%*(-1.0)**K+SM(I,J) 

CONTINUE 
SM(1I,J)=SM(1I,J)+DT(I,J) 

CONTINUE 


CHEEXEKEKEKEEKA KEKE KEKE KEKE RE RETAKE EAE REA TAKER ETE ERE KARE KEE AKT ES 


Cx 
Cx 
Cx 
Cx 
Cx 
Cx 
Cx 
Cx 
Cx 
Ct 
Cx 
Cx 


STT(I,1) - RADIUS OF INTEREST 
STT(I,2) - DEFLECTION 

STT(I,3) - DW/DR 

STT(I,4) - RADIAL MOMENT 

STT(I,5) - SHEAR 

STT(I,6) - TANGENTIAL MOMENT 
STT(C1,7))-) DW 2/0 2 

STT(I,8) - DW*3/DR°*3 

STT(I,9) - RADIAL PLASTIC MOMENT 
STT(I,10)- TANGENTIAL PLASTIC MOMENT 


CEEKKKKKEKEKEKEKEKEK ERE KEKE KEKE REE KEE EEA EKER ERA 


Cx 
C 


Cx 
Cx 


STT(1I,1)=R0 
STT(I,2)=SM(1I,1)/RO 
STT(I,3)=1.0/RO*#¥(SM(1I,2)-STT(1I,2)) 
STT(I,7)=1.0/RO*#(SM(I,3)-2.0*STT(1I,3) ) 
STT(I,8)=1.0/RO*(SM(1I,4)-3.0*STT(I,7)) 
STT(1,4)=-D*(STT(I,7)+NU/ROXSTT(I,3)) 
STT(I,5)=D*¥(STT(1I,8)+1.0/RO*STT(I,7)-1.0/RO*¥2.0% 
STT(1I,3)) 
STT(1,6)=-D*¥(1.0/RO*STT(I,3)+NU*STT(I,7) ) 


THIS IS A STEP TO PREVENT THE INTERIOR POINTS FROM 


Cx HAVING THE EPSILON THROWN IN. 


Ct 


IF (I .EQ. 1) THEN 
RO=B+DEL 

ELSE 
RO=RO+DEL 

ENDIF 
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C 
200 CONTINUE 


C 
RETURN 
END 
C 
Cx 
SSS SSPPSPPSSSSSSESSSSPSSSSSSSSSSSSSSSSOSSSSSLLSPSSSPSSSPSCSCSSL PSP LSS SSS SSS SS - 
Cx 
Cx SUBROUTINE YLD 
Cx 
Cx THIS SUBROUTINE TAKES THE PLATE BENDING SOLUTION FOR THE INPUT 


Cx LOAD INCREMENT AND CALCULATES THE LOAD AT WHICH YIELDING OF THE 
C* OUTER FIBER WILL OCCUR. THE LOAD MATRIX IS THEN INCREASED TO 

Cx THIS VALUE. THIS IS ESSENTIALLY A TIME SAVER, TO AVOID STEPPING 
C* UP THE LOAD IN THE ELASTIC RANGE. 


Cx 
OP PSP PP LSPS LPL LSPS SES 222222925922 929222222°2222 2222222292222 2 2 ff ® & ; 
Cx 
C 
SUBROUTINE YLD( KK, IMAX, LDMAT] ) 
C 
IMPLICIT REAL*8 (A-H,0-Z) 
REAL*8 LDMAT1(50),STRMAX,STTMAX 
INTEGER KK, IMAX, IRMAX, ITMAX 
C 
SINCLUDE: ’COMMON. FOR’ 
C 
Cx 
Cx FIND THE LOCATION ACROSS THE PLATE WHERE THE BENDING MOMENTS 
Cx ARE THE HIGHEST 
Cx 
C 
STRMAX=0.0 
STTMAX=0.0 


DO 122 I=1,NINT(ADJ(2))+1 
IF (ABS(STT(I,4)) .GT. STRMAX) THEN 
STRMAX=ABS(STT(1I,4)) 
IRMAX=I 
ENDIF 
IF CAES(ST2I<( 1,6) ) {~GI-= STIMAX) THEN 
STTMAX=ABS(STT(1,6)) 


ITMAX=I1 

ENDIF 

IF (STRMAX .GT. STTMAX) THEN 
KK=1] 
IMAX=IRMAX 

Evo 
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r22 


Cx 
Cx 
Cx 
Ct 


126 


130 


C 
Cx 


KK=2 
IMAX=ITMAX 
ENDIF 
CONTINUE 


CALCULATE THE EQUIVALENT STRESS AT THE LOCATION WHERE STRESSES 


ARE HIGHEST. COMPARE IT TO THE YIELD STRESS. THEN RATIO UP THE 
LOAD MATRIX ACCORDINGLY 


ERMAX=STT(IMAX,7)*GEOM(5)/2.0 
BTMAX=STT(IMAX,3) *GEOM(5)/2.0/STT(IMAX,1) 
SRMAX=GEOM(6) /GEOM(7)/(1.0-GEOM(3)**2.0)*(ERMAX+GEOM(3) *ETMAX) 
STMAX=GEOM(6)/GEOM(7)/(1.0-GEOM(3)**2.0) *(ETMAX+GEOM(3) *ERMAX) 
SEMAX=SQRT(SRMAX**2.0-SRMAX*STMAX+STMAX*#42. 0) 
R=GEOM(6)/SEMAX 
DO 126 J=1,NINT(ADJ(2)) 
LDMAT1(J)=R*LDMAT1 (J) 
CONTINUE 
WRITE(*,130) LDMAT1(1) 
FORMAT(2X,’LDMAT1 = °,F20.10) 


RETURN 
END 


CHEKEKEEKKKEKKE KEKE KER EK EKA ETEK KEKE REE KEKE KKK KEKE REE EKER KEKE KEKE 


Cx 
Cx 
Cx 
Cx 
Cx 
Cx 
Ct 


SUBROUTINE PLAST 


THIS SUBROUTINE CALCULATES PLASTIC MOMENTS USING A TRAPAZOIDAL 
RULE INTEGRATION SCHEME. THE SUBROUTINE WILL ONLY HANDLE 
LINEAR ELASTIC-PERFECTLY PLASTIC MATERIALS. 


CHEEK KKEKKKEKEKE KEKE KEKE KEE KEKE ERE KEKE KERR KEKE KER KE KEKE EK K EEE EEE 


Ct 
C 


C 


SUBROUTINE PLAST(I) 


IMPLICIT REAL*8 (A-H,O-Z)C 
REAL*8 MRP,MTP,NU 
INTEGER 12 


S$INCLUDE: ’COMMON. FOR’ 
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Cc 
Cx 


SESE SESSESSSLISLSSSSSPSSSSSSSSSESSSSSSSSSSSSSSSSSSSSSSSSSSSS LES SSS SS: 


Cx 
Cx 
Cx 
Cx 
cx 
Cx 
Cx 
Cx 
Cx 
Cx 
Cx 
Cx 
Cx 
Cx 
Cx 


Cx 


VARIABLES USED IN THIS SECTION ARE DEFINED AS FOLLOWS 


DER = TOTAL RADIAL STRAIN INCREMENT 

DET = TOTAL TANGENTIAL STRAIN INCREMENT 

DERE = ELASTIC RADIAL STRAIN INCREMENT 

DETTE] ELASTIC TANGENTIAL STRAIN INCREMENT 

DEL = DEPTH INCREMENT 

st = RADIAL POSITION STATION NUMBER 

TZ = HALF-THICKNESS STATION NUMBER 

DS1 = RADIAL STRESS INCREMENT 

DS2 = TANGENTIAL STRESS INCREMENT 

STR oe MATRIX HOLDING THE TEMPORARY STRESSES AT THE 
GIVEN LOADING INCREMENT (1=RADIAL; 2=TANGENTIAL) 

SET = TEMPORARY EQUIVALENT STRESS - USED FOR YIELD 
CRITERIA 

SE = EQUIVALENT STRESS FROM PREVIOUS INCREMENT 

DEPT = ARRAY CONTAINING THE TEMPORARY STRAIN INCREMENTS 
AT THE GIVEN LOADING INCREMENT (1=RADIAL; 2= 
TANGENTIAL 

R = CONSTANT REFERRING TO THE FRACTION BEYOND YIELDING 
THAT OCCURRED DURING THE COURSE OF THE LAST LOAD 
INCREMENT 

Z = DEPTH 


CHEEK EEK AKER AKER EKA EAA EEE EEE EET 


cx 
C 


C 


C 
Cx 


STT(1,9)= 
STI¢1,10) 
SSUMR=0.0 
SSUMT=0.0 
DEL=GEOM(5)/2.0/(ADJ(6)-1.0) 
Z=0 

E=GEOM(6)/GEOM(7) 

NU=GEOM(3) 


0.0 
=0.0 


DO 500 IZ=NINT(ADJ(6)),1,-1 


OC FSP PPESPSLPPSSEPLLSPSLSESL SSIS SISSSLPSSSSL&PSLSLL LSPS LSE LSE S ESS 2S SS Se S&S @ ® & © @ @ 5 


Cx 
Cx 
Cx 
Cx 
Cx 
C* 
Cx 
Cx 
C* 


THIS SECTION CALCULATES THE STRAIN AND STRESS INCREMENTS 
CORRESPONDING TO THE GIVEN LOAD INCREMENT. IT THEN ADDS THE 
STRESS INCREMENT TO THE PREVIOUS TOTAL STRESS TO CALCULATE THE 
EQUIVALENT STRESS. THE EQUIVALENT STRESS IS COMPARED TO THE 
UNIAXIAL YIELD STRESS TO SEE IF YIELDING HAS OCCURRED. IF SO, 
THE PLASTIC STRAIN INCREMENT IS SET EQUAL TO THE TOTAL STRAIN 
INCREMENT 


OP PSS SSP SESESPLP PSP PPP SSS SPLIP SSP SSSLLSSSELP LLLP SP PPLSEELEP SLL SSDS LSS SDS SSS S&S ® | 


C* 
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DER=-STT(1I,7)*(IZ-1)*DEL 

DET=-STT(I,3)*(IZ-1) *DEL/ST(I,1) 

DERE=DER-DEPT(I,1Z,1) 

DETE=DET-DEPT(I,I1Z,2) 

DS1=E/(1.0-NU*x*2.0)*(DERE+NUXDETE ) 
DS2=E/(1.0-NU%*2.0)*(DETE+NU*DERE) 

STHT (1212, 1)=STR(1,1Z,1)+DS1 

STHRT Gl. 12,2)=STR(1,12,2)+DS2 

SET Glee e5O0n1(STRI(1+1zZ. ))*22-0-STRI(1,12Z,1)*STRT(1,12,2) 


a5 +STRT(I,1Z,2)**2) 
oO 
Cx 
Cx SEE IF YIELD CRITERIA IS EXCEEDED. IF SO, REZERO THE PLASTIC 


Cx STRAIN INCREMENTS AND THE STRESS INCREMENTS. ALSO, RETURN TO THE 
Cx MAIN PROGRAM IF THE OUTER FIBER HAS NOT BEGUN TO YIELD 


Cx 
C 
IF (SET(I,IZ) .LT. GEOM(6)) THEN 
DEPT(I,1Z,1)=0.0 
DEPT (1 ,12,2)=0.0 
Sr (i. 9) =050 
STTCI, 10) =0.0 
IF (IZ .EQ. NINT(ADJ(6))) THEN 
RETURN 
ENDIF 
ELSE 
SET(I,1Z)=GEOM(6) 
DEPT(I,IZ,1)=DER 
DEPT(I,1Z,2)=DET 
[eetSE (1.12) 2LT.. GEOM(G)) THEN 
R=(SET(1,1Z)-GEOM(6) )/(SET(1,1Z)-SE(I,1Z)) 
DEENCl, 12) Ie RSBERT( Liz, )) 
DEPT(1, 1Z,2)=R*DEPT(1, 1zZ,2) 
STRT(I,1IZ,1)=STR(I,1IZ,1)+DS1*(1-R) 
STRT(I,1Z,2)=STR(I,IZ,2)+DS2*(1-R) 
ENDIF 
ENDIF 
c 
Z=Z+DEL 
c 
500 CONTINUE 
c 
Cx 
CHEXESEKEKKEKKEKEKE EEK SEKREK ETEK EKEKAKAEKEAKEKE AEA EKA KK 
Cx 
Cx PERFORM A TRAPAZOIDAL RULE INTEGRATION OF THE STRAIN OVER THE 


Cx PLATE HALF THICKNESS, AND MULTIPLY BY 2 TO GIVE THE TOTAL RADIAL 
C+ AND TANGENTIAL PLASTIC MOMENT INCREMENTS (STT(1I.9) AND STT(I,10)) 
Cx 
CxexExxexEeEEeEKKKEEKEEKEEKEEEEERE KEKE EKER ERE EERE EE EEE EEK 
Cx 
Cc 

Z=-DEL/2.0 

DO 540 IZ=1,NINT(ADJ(6)-12) 
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C 
C 


oe $222 SSS SPSS SSS SS PSP SSS ISSSSPPLE SPS PS SSSSLLE SS SLES PSL LS SES SSS? S&S SF FS F | 


C 
Cc 
Cc 
Cc 


540 


600 


700 


x 


x 
x 
x 
x 


Cx 


C 
Cc 


OF SSS SSS SSSSFSSCSSESESSESLLSPSESLSSSCLSLSESLSESSCSLL ISLS SE SESS L SSS SILLS SSS SS SSDS & 


C 
C 


Cc 


x 
x 


x 


Z=Z+DEL 


DMPR=((DEPT(1,1Z,1)+DEPT(1,1Z+1,1))+GEOM(3)*(DEPT(1I,1Z, 2) 


SDEPT CI Zt), 2)))sZ/2.0 


DMPT=((DEPT(1,1Z,2)+DEPT(I,IZ+1,2) )+GEOM(3)*(DEPT(I,1Z,1) 


STDERR YO), Gori, ])) 2/2. 0 
SSUMR=SSUMR+DMPR 
SSUMT=SSUMT+DMPT 

CONTINUE 


STT(I,9)=2.0*GEOM(6)/GEOM(7)/(1-GEOM(3)**2.0)*DEL*SSUMR 
STT(1,10)=2.0*GEOM(6) /GEOM(7)/(1-GEOM(3) **2.0)*DEL*SSUMT 


WRITE(*,600) I 

FORMAT(2X,’LAT POS = ’,12) 

WHEE (*, 7O00)7 S1T1(1,9)  STT(1,, 10) 
FORMAI( 2X, MHF = “F20. 172," Mie w- 4 prco.l2) 


RETURN 
END 


SUBROUTINE UPDATE 


& 


THIS SUBROUTINE TAKES THE INCREMENTAL STRESSES, STRAINS, 


DEFLECTIONS, MOMENTS, ETC AND ADDS THEM TO THE PREVIOUS TOTALS 
ONCE CONVERGENCE OF PLASTIC MOMENTS HAS OCCURRED. 


SUBROUTINE UPDATE(STT1, IMAX, OUT, LDMAT1, LL, LDMATO ) 


IMPLICIT REAL*¥8 (A-H,0O-Z) 
REAL*8 STT1(50,2), LDMATO(50),LDMAT1(50),NU,DEL,E 
INTEGER NN,OUT,LL 


SINCLUDE: ’COMMON. FOR’ 


Cc 


905 


NN=NINT(ADJ(2) )+1 
E=GEOM(6)/GEOM(7) 
DEL=GEOM(5)/2.0/(ADJ(6)-1.0) 
NU=GEOM({ 3) 


DO 1000 I=1,NN 
DO 905 J=LL, 10 
Stile) =. 10) tod Ow 
CONTINUE 


i22 





Ct 
Cx 
Ct 
Ct 
OR 
On 


OR 
Ct 
On 
Cx 


$80 


C 


LDMATO(1I)=LDMAT1(1)+LDMATO(I) 


DO 980 IZ=1,NINT(ADJ(6)) 


IF THE YIELD CRITERIA HAS BREN EXCEEDED, THEN ADD IN THE 
PLASTIC STRAIN INCREMENTS AND UPDATE THE TOTAL STRESS TO 
FQUAL THE TEMPORARY STRESS CALCULATED IN THE SUBROUTINE 
PLAST 


IF (SE(I,IZ) .GT. GEOM(6)) THEN 
PEG tH RP( LIZ.) (-DRPI( 1, 12,1) 
EPs, 2 -hP (l,l 2,2) 7 DRkP I (1, 12,2) 
Srn( 1,12, )) =SPRI Cl 12,1) 
Ste (l,iZ,2)=STRTI(1,12Z, 2) 
Z=(GEOM(5)/2.0/(ADJ(6)-1.0))#(1Z-1) 
ELSE 


IF THE YIELD CRITERIA HAS NOT BEEN EXCEEDED, THEN UPDATE THE 
STRESSES AND EQUIVALENT STRESSES USING THE ELASTIC EQUATIONS 


STR(I,1IZ,1)=-E/(1.0-NU#*2.0)%*#(ST(1I,7)*(1IZ-1) *DEL- 


= EP(I,IZ,1)+NU*(ST(1I,3)*(IZ-1) *DEL/ST(I,1))) 
STR(I,1Z,2)=-E/(l. O-NU€42.0)*(ST(1,3)4(IZ-1)4DEL/ST(1I,1)- 
a8 EP(I,1Z,2)+NUS(ST(I,7)#(IZ-1) *DEL)) 
Sel IZ)=SCRT(STROI IZ. )s*2-0-STR(IT,1Z,1)% 
- STR(I,1Z,2)+STR(I,1Z,2) *#2) 
ENDIF 
TreeabioCEkP(1,12.,1)) -Gl. GEOM(S) .OR. ABS(EP(1,1Z,2)) 
Yy .GT. GEOM(9)) THEN 
OUT=1 
ENDIF 
CONTINUE 
1000 CONTINUE 
LL=2 
RETURN 
END 
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Cc 


Cx 
CEEKEEKEEEKEKEKEKEK EK KAKA KEKE EKER EEE KEE EEE EEE EEE 
Cx 
Cx SUBROUTINE RESUL 
Cx 
Cx STORE THE RESULTS AT THE END OF EVERY 10 LOAD INCREMENTS IN 
cz THE FILE OUT.DAT. ALSO SCREEN DUMP THESE RESULTS. 
Cx 
CHEEXKKERKEEE EKER EAA ERA KEKE KEKE AKA EKER TAREE KEKE EKA EKEKEKE 
Cx 
C 
SUBROUTINE RESUL(MM,STT1,COUNT, LDMATO) 
04 
IMPLICIT REAL*8 (A-H,0-Z) 
REAL*8 STT1(50,2),LDMATO(50) 
INTEGER NN 
SINCLUDE: ’ COMMON. FOR’ 
Cc 
NN=NINT(ADJ(2)) +1 
DEL=(GEOM(1)-GEOM(2))/ADJ(2) 
Cc 
IF (COUNT .EQ. 1.0) THEN 
WRITE(15,100) 
WRITE(*,100) 
100 FORMAT (///, ’ ¥¥¥¥EEEREEEKEEEEKEREREKEAKAA AAA EEEE’ , 
+ "EREKEKEKREKER EAA KEKE EAE EEE’ //) 
WRITE(15,120) MM, LDMATO(1)/DEL 
WRITE(*,120) MM, LDMATO(1)/DEL 
120 FORMAT(2X,’ LOAD INCR = ’°,I3,’ LOAD = *,F10.4,7) 
c 
WRITE(15,130) 
WRITE(*,130) 
130 FORMAT(2X, ’RADPOS’,2X, DEFLECT °,2X,’SLOPE °,2X,’RAD MOM ’ 
+ ,2X,’ MRP °,2X,’ TAN MOM '°,2X,’ MTP ’,’ # £=SHEAR oh) 
DO 160 I=1,NN 
STM1=ST(1I,4)+ST(I,9) 
STM2=ST(I,6)+ST(I,10) 
WRITE(15,140) ST(I,1),-ST(I,2),ST(1I,3),STM1,ST(I,9), 
+ SIMZ,STCI,10),ST(1,5) 
WRITE (£, 140) STC1,1),—-ST(1, 2).S57(1,3),STIM1,ST(1,9), 
+ STMZ,ST(1,120),STC1,5) 
140 BORMAT (2X8 5-2,2%, F105, 2X87 -4,2%,F10.0,2X,F5.0,2x%,F10.0, 
+ 2xX,F5.0,2X,F10.0) 
160 CONTINUE 
Cc 
WRITE(15,163) 
WRITE(*,163) 
163 Cy eG 00 Nr i ee ees eee are Re 2 .; 
+ 4 Stier yi oad SMe G eninge aR on) 
DO 1000 I=1,NN 
IF (ST(I,9) .NE. 0.0 .OR. ST(I,10) .NE. 0.0) THEN 
WRITE(15,170) ST(I,1) 
WRITE(*,170) ST(I,1) 
170 FORMAT(/,2X, RADIAL POSITION = ’°,F5.2) 


Pod 








WRITE(15,165) 
WRITE(*, 165) 
165 FORMAT(/,2X,° DEPTH °,2x,’ ERP *.2x,’ SRAD 
+ 2y ETP “2% | STAN OX SE ip) 
DO 200 IZ=1,NINT(ADJ(6) ) 
IF (SE(I,IZ) .GT. GBOM(6)) THEN 
Z=(GEOM(5) /2.0/(ADJ(6)-1.0))*(1IZ-1) 
WRITE(15,180) Z,EP(1,1Z,1) JSTR(I,1Z,1),8P(1,12Z,2), 
+ STR(1L,1Z;2),SE(1,1Z) 
WRITE(*,180) Z,EP(I,1IZ,1),STR(I,1IZ,1),B8P(I,1IZ,2), 
+ STR(1,1Z,2) J/SE(1,1Z) 
180 FORMAT(2X,F6.4,3X,F10.8,3X,F8.1,3X,F10.8,3X,F8.1, 
+ SX.6S.1) 
ENDIF 
200 CONTINUE 
ENDIF 
1000 CONTINUE 
ENDIF 
RETURN 
END 


COMMON GEOM(9),ADJ(7),ST(50,10) 

COMMON /PLSCM/ EP(50,50,2),STT(50,10), DEPT(50,50,2) ,STR(50,50,2), 
7 Si 000,50) ,S5STRT(50,50,2),SET(50,50) 

REAL#8 ADJ,ST,EP,STT,DEPT,STR,SE,STRT,SET 


LS 
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